Celera Assembler

Der Celera Assembler, ein Genom-Assembler, wurde ursprünglich von dem Unternehmen Celera entwickelt und wird nun als Open-Source-Projekt weitergeführt. Er wird dazu genutzt, aus vielen kurzen genomischen Fragmenten, die durch eine Whole-Genome-Shotgun-Sequenzierung gewonnen wurden, eine lange zusammenhängende Sequenz zu rekonstruieren.

Dabei werden mit einem Seed-and-Extend-Ansatz Überlappungen zwischen den bereinigten Fragmenten gesucht. Daraus wird der Overlap-Graph konstruiert. Eine Spanning-Forest-Heuristik setzt die Fragmente zu Unitigs zusammen. Diese wiederum werden vom Greedy-Path-Merging-Algorithmus mit Hilfe der Mate-Pair-Informationen und dem Contig-Mate-Graph zu Scaffolds arrangiert. Die Lücken zwischen den Contigs können eventuell mit bisher ungenutzten Fragmenten und Mate-Pairs heuristisch gefüllt werden. Am Ende wird die Konsensussequenz durch ein Multialignment aller Fragmente über den gefundenen Scaffolds berechnet.

Theoretische Betrachtungen

Ein Assembly lässt sich durch folgendes Problem beschreiben.

Es sei eine Menge $F = {f_2, f_2,\dotsb,f_n}$ von Strings gegeben. Gesucht ist ein String S mit folgenden Eigenschaften:

1. Jedes f_i ist Substring von S, $\forall f_i \in F : f_i \in S$
2. S ist minimal, $\forall S'$ mit 1. gilt $|S| \le |S'|$

Dieses Problem ist als Shortest Common Superstring bekannt. Es ist NP-vollständig.

Im biologischen Kontext kommt verschärfend hinzu, dass die Sequenzen fehlerhaft sein können und stets zwei mögliche Orientierungen eines Fragmentes zu berücksichtigen sind. Außerdem ist es aus biologischer Sicht nicht immer sinnvoll eine minimale Sequenz zu ermitteln, da Fragmente aus repeathaltigen Regionen zu überkomprimierten Ergebnissen führen können. Das und die kombinatorische Komplexitiät bewirken, dass das hier vorgestellte Verfahren stark heuristisch ist.

Verfahren

Der Celera Assembler implementiert den Whole-Genome-Shotgun-Ansatz. Der Assembler lässt sich nach [MSD+00] in mehrere, aufeinander aufbauende Stufen (Abbildung 2) einteilen, welche alle heuristisch sind. Das Herzstück des Assemblers bildet der Greedy-Path-Merging-Algorithmus.

Eingabedaten

Die Gewinnung der Sequenzfragmente geschieht mit Hilfe der "Double Barrel" Shotgun-Sequenzierung. Dabei wird ein Klon von beiden Seiten ansequenziert und man erhält ein sogenanntes Mate-Pair von Reads. Dabei verwendete Klone haben typische Längen von 2 kb, 10 kb, 50 kb und 150 kb und die sequenzierten Reads eine Durchschnittslänge von 550 Basen. Die ursprünglichen Reads werden daraufhin auf ein Intervall mit Basen 98-prozentiger Sicherheit gekürzt. Die den gelesenen Basen zugeordneten Qualitätswerte werden im weiteren Prozeß nicht berücksichtigt. In den Reads enthaltene Teile von Klonierungs- oder Sequenzierungsvektoren werden ebenfalls radikal entfernt. Die verbleibenden Stücke hoher Qualität sind die Eingabe des Assemblers. Diese werden im folgenden als Fragmente oder Reads bezeichnet.

Screener

Die Sequenzfragmente werden zunächst nach bekannten Repeats durchsucht. Fragmente, welche Repeatmuster enthalten, werden maskiert und später gesondert betrachtet. Die bereinigten Fragmente bilden die Eingabe für den nächsten Schritt.

Overlapper

Der Overlapper sucht Überlappungen zwischen Fragmenten. Eingabe dieser Phase sind die beim Screening bereinigten Fragmente. Ausgegeben werden alle signifikanten Überlappungen zwischen den Reads. Für jedes Fragment f_i ist also zu bestimmen, wie es mit irgendeinem anderen Fragment f_j oder dessen reversen Komplement $\overline{f_j}$ überlappt. Zwei beliebige Fragmente f_i und f_j können auf verschiedene Weise überlappen (Abbildung 4).

Abbildung 4: Mögliche Überlappungen von Fragmenten

Das heißt, sie überlappen an den Enden, eines ist im anderen enthalten oder sie überlappen nicht. Alignment durch dynamische Programmierung ist bei dieser großen Anzahl Fragmente zu langsam. Beim humanen Genom wären rund 27 Mio. Reads paarweise zu vergleichen. Hier sind außerdem nur Überlappungen mit einem hohen Alignmentscore, das heißt mit mehr als 96 % Identität, relevant. Daher wird ein BLAST-ähnlicher "Seed-and-Extend"-Ansatz verfolgt (Abbildung 5).

Abbildung 5: Suche nach Seeds

Alle Fragmente werden zunächst konkateniert, sei also $H = f_1f_2 \dotsb f_n$. Anschließend sucht man für alle Fragmente f_i exakte Matches der k-mere von f_i auf der Sequenz H, die sogenannten Seeds. Der Parameter k wird hier zwischen 18 und 22 gewählt. Nun versucht man die Seeds durch Banded Alignment zu erweitern. Die Laufzeit ist linear und nur wirklich gute Überlappungen zwischen f_i und dem entsprechenden f_j werden gefunden (Abbildung 6).

Abbildung 6: Erweiterung eines Seeds durch Banded Alignment

Wurde eine Überlappung zwischen fi und fj gefunden, kann das drei Ursachen haben. Die Fragmente überlappen tatsächlich auf der Originalsequenz wie f1 und f2 in Abbildung 7, die Enden der Fragmente stammen aus sich wiederholender Sequenz (Repeat) wie f3 und f4 in Abbildung 7 oder die Sequenzgleichheit ist rein zufällig. Letzteres wird dadurch vermieden, dass nur Überlappungen einer Länge von mindestens 40 bp betrachtet werden. Die durch Repeats erzeugten Overlaps können vermieden werden, indem man vorher die Fragmente penibel nach bekannten Repeats screent und maskiert oder indem man k fest wählt, für jedes k-mer die Häufigkeit bestimmt und extrem oft auftretende k-mere ignoriert. Die berechneten Überlappungen zwischen den Fragmenten werden nun in einem Graphen repräsentiert.

Abbildung 7: Wahre und falsche Overlaps

Der Overlap Graph wird wie folgt definiert: für jedes Fragment fi gibt es zwei Knoten, einen Startknoten si, einen Endknoten ei und eine gerichtete Kante (si, ei), um die Orientierung des Reads darzustellen. Jede Überlappung wird durch eine ungerichtete Kante dargestellt. Die Overlap-Kanten inzidieren mit Knoten an den Enden, an denen die entsprechenden Fragmente überlappen (Abbildung 8).

Abbildung 8: Overlap und Overlap-Kante

Die Länge einer Kante ist die negative Länge der Überlappung. Das heißt, große Überlappungen bedeuten kurze Kanten. Im folgenden Beispiel sei F = {f1, f2, . . . , f7} die Menge der Fragmente mit Überlappungen wie in Abbildung 9.

Abbildung 9: Overlaps im Beispiel

Hier überlappen die ersten 40 Basen von Read f1 mit den ersten 40 Basen von f4, die letzten 300 Basen von f1 mit den ersten von f6 und so fort. Mit der Definition ergibt sich daraus der Overlap Graph in Abbildung 10.

Unitigger

Der Unitigger nutzt den Overlap Graph, um die Fragmente zunächst anzuordnen, das heißt allen Knoten Koordinaten zuzuweisen (Layout). Als Eingabe dienen die Überlappungen zwischen Fragmenten aus der Overlap-Phase. Eine einfache, auch hier verwendete Heuristik ist die des spannenden Waldes. Zunächst werden alle Read-Kanten markiert. Die Overlap-Kanten sortiert man aufsteigend nach Länge und fügt solange die kürzesten Kanten hinzu, bis in jeder Zusammenhangskomponente ein spannender Baum entstanden ist. Overlap-Kanten, die einen Kreis schließen würden, werden ignoriert. Im Graph aus dem vorigen Beispiel beginnt man mit der Kante (e3, e4), die mit -420 die "kürzeste" ist. Die Kanten der Längen -300, -280, -260 und -200 kommen dazu. Die Kante mit l = -180 würde einen Kreis schließen. Mit der Kante der Länge -140 ist der Spannbaum fertig (Abbildung 11).

Aus dieser Teilmenge der Kanten ergibt sich eine relative Anordnung der Reads dieser Zusammenhangskomponente. An der Achse in Abbildung 12 stehen die, den Start- und Endknoten zugeordneten, Koordinaten. Solch ein vermeintliches Multialignment heißt Contig.

Durch Repeats erzeugte falsche Overlap-Kanten können zum falschen Zusammensetzen eines Contig (Misassembly) führen. Aus folgender Situation ergibt sich dieser Overlap-Graph.

Ein Spannbaum, der die falsche Kante e enthält, würde eine falsche Anordnung der Reads erzeugen. Lässt man die falschen Kanten außer Acht ergibt sich natürlich die richtige Anordnung. Es gibt aber in diesem Fall kein Layout, das mit allen gefundenen Überlappungskanten konsistent ist. Daher definiert man: Ein Unitig (unique contig) ist ein Contig, das mit allen Kanten des Overlap Graphen konsistent ist. Reads, die nicht zu einem Unitig gehören, werden im Folgenden nicht verwendet. Nun kann es vorkommen, dass ein längeres Stück Sequenz wiederholt auftritt. Das führt dazu, dass ein Contig mit allen Überlappungskanten des Overlap Graphen konsistent ist, dessen Fragmente aber nicht aus einer einzigartigen Region der Ursequenz stammen. Man definiert weiterhin: Ein U-Unitig (unique unitig) ist ein Unitig, das einzigartig auf der Ursequenz ist, also nicht teilweise oder vollständig in einer sich wiederholenden Region liegt. U-Unitigs werden mit Hilfe statistischer Betrachtungen identifiziert. Es wird angenommen, die Ankunft der Fragmente wäre poissonverteilt. Sei R die Anzahl Reads, G die Länge der Ursequenz und u die Länge eines Unitigs. Dann ist R/G die durchschnittliche Anzahl Reads pro Base und $\lambda = \frac{u \cdot R}{G}$ die erwartete Anzahl Reads pro Unitig (Erwartungswert). Die Wahrscheinlichkeit, dass ein Unitig k Fragmente enthält, ist mit der Poisson-Verteilung $\frac{\lambda^{k}} {k!} e^{-\lambda}$ (linke Kurve in Abbildung 13). Entstand das Unitig aus den Reads zweier Repeats, verdoppelt sich der Erwartungswert und die Wahrscheinlichkeit ist $\frac{2\lambda^{k}} {k!} e^{-(2\lambda)}$ (rechte Kurve in Abbildung 13).

Literatur

• "Algorithms in Bioinformatics II", Kapitel 17, Prof. Dr. Daniel Huson, Universität Tübingen, 2004 (engl.) (PDF-Datei; 1,49 MB)
• Der Celera Assembler", Seminararbeit, 2005 (PDF-Datei; 890KB)

Weblinks

• Projektseite auf SourceForge