Leyendo el diagrama de Ramachandran

Hola,
la semana pasada les planteé a mis alumnos de la Licenciatura en Ciencias Genómicas el problema de asignar estructura secundaria a los residuos de una proteína, una vez medidos los ángulos dihedros del esqueleto peptídico psi y phi. Para esta tarea lo lógico es usar el diagrama de Ramachandran, como éste tomado de la wikipedia:

La solución naive consiste en delimitar las regiones del diagrama por medio de condiciones tales como 'si phi > X  AND psi < Y' entonces un residuo está en estado Z. Sin embargo, para tener mucha precisión sería necesario tener muchas condiciones como ésta, dada las formas irregulares del diagrama.


Una solución interesante, propuesta por F.Peñaloza, consiste en realmente leer el mapa, como lo hacemos nosotros. Lo que él hizo fue obtener un mapa de Ramachadran cuadrado con unos cuantos colores conocidos,

Diamagra de Ramachandran RGB de F.Peñaloza.

Diagrama anterior con los colores RGB empleados

para luego convertirlo a texto ASCII (con algo como text-image). Con este mapa en memoria, es sencillo obtener la estructura secundaria para unas coordenadas phi,psi dadas. 

Aprovechando el módulo GD de perl podemos leer la imagen directamente, como muestra el siguiente código:

 use strict;  
 use GD;  
   
 my $RAMAPLOT = 'ramachandran.png';  
 my $plot = GD::Image->newFromPng($RAMAPLOT);  
 my ($maxpsi,$maxphi) = checkPlot($plot);  
   
 my $inputPSI = 160;   
 my $inputPHI = -45;   
 printf("El estado de estructura secundaria para %s,%s es: %s\n",  
    $inputPSI,$inputPHI, getSimpleSS($maxpsi,$maxphi, $inputPSI,$inputPHI));  
   
 sub getSimpleSS  
 {  
   my ($maxpsi,$maxphi,$psi,$phi) = @_;  
   my %colorSS = (  
     '236,30,36','E',  
     '164,74,164','E',  
     '156,218,236','H',  
     '4,162,236','H',  
     #'180,230,28','L'   
     );  
   
   if($psi > 180 || $psi < -180 || $phi > 180 || $phi < -180)  
   {  
     die "# valid psi/phi are [-180,+180]\n";  
   }  
   $psi = (0.5 - (0.5 * $psi/180)) * $maxpsi;  
   $phi = (0.5 + (0.5 * $phi/180)) * $maxphi;  
     
   my @RGB = $plot->rgb( $plot->getPixel($phi,$psi) );  
   return $colorSS{"$RGB[0],$RGB[1],$RGB[2]"} || 'C';  
 }  
   
 sub checkPlot  
 {  
   my ($plot,$print_pixels) = @_;  
   my ($maxphi,$maxpsi) = $plot->getBounds();  
   if(!$print_pixels){ return ($maxphi,$maxpsi) }  
     
   foreach my $phi ( 0 .. $maxphi )   
   {  
     foreach my $psi ( 0 .. $maxpsi )   
    {  
       my @RGB = $plot->rgb( $plot->getPixel($phi,$psi));  
       print "$RGB[0] $RGB[1] $RGB[2]\n";  
    }  
   }  
   return ($maxphi,$maxpsi);  
 }  

Un saludo, Bruno

HMMER 3.0, HMMER 2.3.2 or PfamScan?

Last time that I annotated Pfam domains into footprintDB database I used the program hmmpfam from the HMMER 2.3.2 software package. But now, many things have changed, Pfam version has moved from 23 to 26, and the current HMM file can't be used directly with HMMER 2, it needs to be converted with a simple tool of HMMER 3 (hmmconvert -2 Pfam-A.hmm > Pfam_ls_26).

HMMER 3 is a nice software tool that is hundreds of times faster than its predecesor, it takes 20 minutes in my Quad-Core computer the same calculation that took like 2 hours in a 28 node cluster.

So I have decided to move to modern times, but cautiously, because last time I tried HMMER 3 I had not wanted results, so I have done my own benchmark that I'm going to explain...

First, I downloaded a test set of protein sequences from 3Dfootprint, as I work with DNA binding proteins, I downloaded all of them from this archive (currently 2007 FASTA sequences).

To calculate the Pfam domains I used the last version of HMMER 3.0 from http://hmmer.janelia.org/software, my old version of HMMER 2.3.2 (something similar can be found here: http://hmmer.janelia.org/software/archive) and pfam_scan.pl script used by Pfam team to create their database in the Sanger Institute. Also I downloaded the last version of Hidden Markov Models from Pfam version 26 and converted it to use with HMMER 3 (hmmpress Pfam-A.hmm) and with HMMER 2 (hmmconvert -2 Pfam-A.hmm > Pfam_ls_26).

Then I started the testing with the 3 programs, with and without using thresholds in the HMMER param options:
HMMER 2 with thresholds: hmmpfam --acc --cut_ga Pfam-A.hmm protein_sequence_complexes.faa > protein_sequence_complexes.hmmscan
HMMER 2 default: hmmpfam --acc Pfam-A.hmm protein_sequence_complexes.faa > protein_sequence_complexes.hmmscan
HMMER 3 with thresholds:  hmmscan --acc --cpu 8 --notextw --cut_ga -o protein_sequence_complexes.hmmscan Pfam-A.hmm protein_sequence_complexes.faa
HMMER 3 default: hmmscan --acc --cpu 8 --notextw -o protein_sequence_complexes.hmmscan Pfam-A.hmm protein_sequence_complexes.faa
PfamScan (default thresholds): pfam_scan.pl -align -cpu 8  -hmm Pfam-A.hmm -fasta protein_sequence_complexes.faa -outfile protein_sequence_complexes.pfamscan

The final conclusion is that I'll use HMMER 3 with thresholds, it's because the calculation time is 200 times faster that HMMER 2 (Figure 1) and both retrive more or less the same number of domains for the main transcription factor families (Figure 2).

It's remarkable that HMMER 3 without thresholds is very much sensitive, detecting near the double number of domains than the rest of the techniques (Figure 1), but most of them are undesired domains that interfere with the identification of the important ones.

HMMER 2 results with and without thresholds are comparable (Figure 1), both of them detect most of the transcription factor domains (Figure 2), even a bit more than HMMER 3 with thresholds, without including spurious domains even without thresholds (Figure 1).

PfamScan detects less domains than the rest of the techniques (Figure 1), although it uses HMMER 3 internally, this is because it doesn't annotate overlapping domains, but also because it has very strict thresholds that in many cases fail to detect real transcription factor domains (Figure 2). We have noticed this problem in a particula recent study of the transcription factor YY1 (1UBD chain C), if we search in Pfam webserver its sequence (chain C) we obtain only 3 of the 4 real Zinc Finger domains, we must find the 4th Zinc Finger domain among the 'insignificant Pfam-A Matches'.

I hope these results help to decide to people like me dubbing among moving to HMMER 3, use PfamScan or continue using HMMER 2.

Figure 1. Statistics of several parameters with the 5 calculation methods.

Figure 2. Numer of retrieved domains for different transcription factor families with the different methods.

BIFI 2012 - V Congreso Internacional - Dianas proteicas: Descubrimiento de Compuestos Bioactivos

Del 1 al 4 de febrero de 2012 se celebrará en Zaragoza la 5º Congreso Internacional del Instituto de Biocomputación y Física de Sistemas Complejos (BIFI).

Este año el tema central será el Descubrimiento de Fármacos, cubriendo desde los pasos iniciales de investigación en laboratorio hasta los estudios preclínicos: nuevas dianas proteicas, validación de dianas, nuevas metodologías y herramientas de caracterización tanto estructural como funcional y cribado computacional de moléculas. La conferencia servirá de punto de encuentro de investigadores del campo del descubrimiento de fármacos, donde se discutirán los avances más recientes y retos futuros.

Nuestro laboratorio asistirá al evento con un seminario titulado "Protein-DNA interface prediction techniques: performance and potential in protein engineering" y un póster titulado: "In vivo DNA binding pattern of Rex-1 in mouse embryonic stem cells" realizado en colaboración con el Departamento de Veterinaria de la Universidad de Zaragoza.

English version:

The V International Conference of the Institute for Biocomputation and Physics of Complex Systems (BIFI) on February 1-4, 2012.

The meeting will be an international conference on Drug Discovery from a protein perspective, covering most of the initial steps in drug discovery and preclinical studies (new protein targets, protein target validation, new methodologies and tools for structural and functional characterization, experimental and computational high-throughput screening, etc.). We wish the conference to represent a venue for gathering active researchers on drug discovery, with strong roots in the scientific and academic communities to discuss recent developments and future challenges in the field.

Our laboratory will participate in the event with a talk titled "Protein-DNA interface prediction techniques: performance and potential in protein engineering" and a poster titled: "In vivo DNA binding pattern of Rex-1 in mouse embryonic stem cells" in collaboration with the Veterinary Department of the University of Zaragoza.

Apuntes sobre las XI Jornadas Bioinformáticas

Hola,
tras el anuncio de hace unas semanas, voy a comentar un poco mis impresiones sobre las XI Jornadas de Bioinformática. Después de la sesión inaugural del lunes, a su vez precedida por el simposio de estudiantes, a las que no pude asistir, el día 24 realmente empezaron mis jornadas. Como valoración general, creo que los asistentes hemos tenido ocasión de aprender y de discutir sobre los problemas del campo, y hemos tenido la oportunidad de escuchar charlas muy buenas. Para los que no pudieron venir, ahí va mi resumen de las charlas a las que asistí y notas sobre algunos pósters.

Modelo 3D de la region cromosómica ENCODE ENm008, que contiene el locus de la α-globina, adaptado de www.ncbi.nlm.nih.gov/pmc/articles/PMC3056208.

24 de Enero
Hoy hemos podido escuchar un montón de charlas sobre temas variopintos, y me ha llamado la atención la abrumadura presencia de la palabra mágica NGS (Next Generation Sequencing) en muchas de ellas. También en los pósters expuestos ha habido muchos ejemplos de la aplicación de estas herramientas, sobre todo de RNAseq.

De lo que he visto muy poco hoy ha sido de ChIPseq, tan sólo el póster de Ionas Erb con el software Pro-Coffee para alinear secuencias de DNA de promotores, publicado en NAR.

Otro póster que me llamó la atención fue el trabajo de Minoche sobre la evaluación sistemática de los errores típicos de la plataforma de secuenciación Solexa, publicada en Genome Biology.

Sonia Tarazona me explicó su póster sobre RNAseq y me invitó a probar el  software Qualimap, que permite evaluar el efecto del coverage sobre la interpretación de las diferencias de expresión medidas en RNAseq.

Leo Mirny nos explicó en detalle la técnica de Hi-C para localizar regiones de cromatina cercanas en el espacio celular y cómo en su grupo han usado este tipo de datos para entender el empaquetamiento del núcleo de levaduras y de Homo sapiens, construyendo una matriz que se parece mucho a un mapa de contactos de proteínas, algo que también por la tarde explicó Davide Bau en su trabajo sobre la transcripción y la estructura de la cromatina en un locus de alfa globina.

Juan Ramón González nos comentó los métodos predominantes actuales para la normalización de conteos de lecturas RNAseq (TMM, EDAseq y CQN) y presentó los problemas que tienen las distribuciones de Poisson (con un parámetro libre) y la binomial negativa (con dos) para modelar algunos datos reales, y mostró ejemplos convincentes del uso de la de Poisson-Tweedie como lo mejor de los dos mundos, con un tercer parámetro para elegir según el caso el mejor modelo estadístico, dada la dispersión de los datos reales de esta teconología. Propone su paquete de Bioconductor/R tweeDEseq como herramienta para esta tarea.

Eva María Novoa nos dió una clase magistral sobre el uso de codones en procariotas y eucariotas, presentando evidencia de la importancia de las enzimas UMS (en proka) y hetADATS (adenosine deaminasas de euka), que modifican terceras bases de los tRNAs, para explicar las diferentes frecuencias de codones en todos los bichos conocidos. Su trabajo se publica en Cell.

Nacho Medina nos deslumbró con la capacidad de su equipo del CIPF para crear una tubería de análisis de datos de NGS que nos permitirá como usuarios hacerlo todo en sus servidores en "tiempos de minutos", aprovechando la optimización que han hecho de los distintos algoritmos y del hardware subyacente, que incluye, si no recuerdo mal, CPUs, GPUs dedicadas y discos de estado sólido. Me llamó mucho la atención su navegador genómico HTML5, que se puede probar en http://genomemaps.org.

Tomas Marques nos volvió a hablar de primates, ya lo había hecho en Málaga en el 2010. Esta vez trató de convencernos de la importancia de mirar con lupa los datos de NGS, en su caso de ensamblaje genómico, de usar el software adecuado para nuestros objetivos, y de hacer el control de calidad en casa, sin delegarlo. Menciona que en su labo usan GATK como software para variant calling tras haberlo comparado con otros.

Tanya Vavouri nos explicó, si lo entendí bien, que los espermatozoides humanos maduros conservan sólo un 4% de los nucleosomas, pero que justo esos pueden ser muy importantes para pasar información epigenéticas al nuevo cigoto, porque se correlacionan con picos de %GC muy cercanos a promotores.

Javier Macia nos mostró ejemplos de cómo modelar circuitos electrónicos a base de puertas lógicas implementadas con células de levadura modificadas.

Ya hacia el final del día Toni Giorgino nos mostró dinámicas moleculares espectaculares de un dominio SH2 y su ligando, y Pablo Minguez publicó un montón de resultados sobre la conservación de sitios en proteínas eucariotas que pueden sufrir modificaciones postraduccionales. Lo siento, a las dos últimas charlas no me pude quedar.

25 de Enero
El tercer y último día del congreso arrancó con una conferencia plenaria de Luis Serrano donde nos resumió los estudios de su grupo sobre Mycoplasma pneumoniae, un parásito bacteriano con genoma extremadamente reducido que se puede cultivar en el laboratorio. Su charla se puede resumir  como la aplicación de todas las herramientas bioinformática, genómicas, proteómicas y metabolómicas disponibles para tratar de caracterizar la biología de este bicho, que tiene solamente unos 10 factores de transcripción homólogos de otros conocidos en otras especies. Por destacar algo de una charla muy densa pero amena, Luis habló que encontraban una correlación <0.50 entre los niveles de expresión génica y las cantidades de proteína detectadas por MS en distintas condiciones. Esta observación no es muy novedosa, pero sí la explicación que proponía, basada en la divergencia de las secuencias Shine-Dalgarno en los mensajeros, que motivarían menores afinidades de los ribosomas y por tanto menores tasas de traducción.

Antonio Mérida presentó el software Sma3s para la anotación de genomas, que comparó con otros programas como Blast2GO. Roderic Guigó apuntaba  tras la charla que la anotación de un genoma actualmente debe incluir no sólo genes codificantes sino también RNAs reguladores, por ejemplo.

Paolo Ribeca dió la primera charla del día dedicada al tema estrella (NGS). En esta conferencia Paolo repasó los principales problemas del mapeo de lecturas (reads) sobre un genoma de referencia y las limitaciones de los principales programas (Bowtie,BWA,SOAP,MrFAST,MrsFAST) a la hora de hacer búsquedas exhaustivas y flexibles sobre posibles dianas genómicas, algo que su propia plataforma GEM parece haber resuelto y acelerado considerablemente. Su mensaje de precaución es que el usuario de este tipo de software debe conocer con cierta precisión cómo funciona el programa que va a usar y qué limitaciones tiene, en vez de confiar en el programa a ciegas y dejar que tome decisiones, no siempre transparentes, por ti.

Darío Guerrero mostró datos sobre la validación de una tubería de preprocesamiento de lecturas NGS y el ensamblaje de datos de RNAseq, con seqtrimnext y fulllengther, respectivamente. Uno de los programas con los que comparó sus resultados fue Mira.

Beatriz García nos contó el desarrollo de software de aprendizaje automático para la asignación de secuencia de proteínas sin anotar a rutas metabólicas.

Ya en la tarde, Patrick Aloy nos contó un proyecto de su grupo que reconstruye una red de interacciones de proteínas implicadas en la enfermedad de Alzheimer y su uso, junto con una base de datos de fármacos y sus efectos terapéuticos, para predecir el efecto de nuevos compuestos para su tratamiento, así como para replantearse el uso de otros. Parte de este trabajo está publicado aquí.

Ana Rojas nos explicó como su grupo había reconstruido la filogenia de la superfamilia de proteínas RAS, encontrando por el camino los residuos funcionales responsables de las diferencias funcionales de las diferentes familias.

Mar Gonzàlez  nos resumió resultados de su reciente artículo sobre la variabilidad del splicing alternativo en poblaciones humanas.

Alberto Pascual García nos mostró cómo, a partir de datos de presencia/ausencia de rRNA en muestras de diferentes ambientes, se puede inferir la composición bacteriana en un ambiente y por medio de una aproximación basada en la arquitectura de las redes resultantes estudiar si los diferentes géneros tienden a agregarse o segregarse.

Hernán Dopazo nos contó resultados de un trabajo suyo que está en revisión donde sostiene que la composición de los genomas, en cuanto a elementos como genes, rRNAs, promotores, elementos repetidos, etc, se distribuye de manera parecida a la distribución de especies que se observa en ecosistemas naturales, en un proceso donde aparentemente la selección tiene poco que decir.

En la última charla que pude escuchar Jaime Huerta nos contó los progresos de su grupo para hacer filogenias anidadas, que se pueden entender como un proceso recursivo donde vamos refinando el árbol inicial recalculando de manera recursiva la topología de las ramas a medida que vamos de la raíz a las hojas, añadiendo nuevos genes ortólogos a medida que avanza el proceso.

Entre los pósters del segundo día tomé nota del servidor iLOOP para la predicción de interacciones proteína-proteína, del software TAPyR para el alineamiento de reads largos como los de 454 y del programa Pyicos para el procesamiento de datos de ChIPseq.

PD Una oferta de trabajo que se publicó en el congreso:

The Evolutionary Genomics Group in the Comparative and Computational Genomics program of the IBE (http://www.ibe.upf-csic.es/) is willing to recruit a PhD student. More information is available in the attached document. For queries, please contact Tomás Marquès-Bonet (tomas.marques@upf.edu)

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$