Post on 22-Jan-2021
transcript
MASTERARBEIT
Multiples Sequenzalignment mit HiddenMarkov Modellen
durchgefuhrt amStudiengang Informationstechnik und System–Management
an derFachhochschule Salzburg
vorgelegt von:DI(FH) Roland J. Graf
Studiengangsleiter: FH-Prof. DI Dr. Gerhard JochtlBetreuer: Univ. Doz. Dr. Stefan Wegenkittl
Salzburg, Oktober 2010
Eidesstattliche Erklarung
Hiermit versichere ich, Roland J. Graf, geboren am 6. Mai 1964, dass die vorliegendeMasterarbeit von mir selbstandig verfasst wurde. Zur Erstellung wurden von mir keineanderen als die angegebenen Hilfsmittel verwendet.
0910581005Roland Graf Matrikelnummer
ii
Informationen
Vor- und Zuname: DI(FH) Roland J. GrafInstitution: Fachhochschule Salzburg GmbHStudiengang: Informationstechnik & System-ManagementTitel der Masterthesis: Multiples Sequenzalignment mit Hidden
Markov ModellenBetreuer an der FH: Univ. Doz. Dr. Stefan Wegenkittl
Schlagworter
1. Schlagwort: Multiple Sequence Alignment2. Schlagwort: Hidden Markov Model3. Schlagwort: Alignment Scoring
Abstract
This thesis describes the development of a method for performing Multiple SequenceAlignments (MSAs) using Hidden Markov Models (HMMs).The developed solution implements a dynamic, iterative process for generating theMSA. After training the model using pre-aligned sequences (profile-MSA), sequencesare ranked and aligned with the model individually. The highest ranked sequences arethen merged with the profile-MSA, before an alternative algorithm (ClustalW) is usedto make partial improvements to the alignment. In a further step, the generated MSAis prepared such that it allows for an iterated profile-extension scheme much alike PSI-BLAST where new HMMs are generated on the basic of top scoring hits and theiralignments to the original MSA.The quality of alignments and the quality of progressive growing MSA is evaluated bygraphical and numerical comparisons. Both probability density functions of Z-scores aswell as the scattering of Z-scores are used to evaluate the discriminatory power of themodel. Finally shows the change of entropy, whether additional sequences increase theaverage information content of the emission matrix, and whether these changes affectthe results.
iii
Inhaltsverzeichnis
Eidesstattliche Erklarung ii
Informationen iii
Schlagworter iii
Abstract iii
Abbildungsverzeichnis v
Tabellenverzeichnis vi
1 Einfuhrung 1
1.1 Problemstellung und Motivation . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Uberblick . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Allgemeine Grundlagen 5
2.1 Proteinstruktur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Sequenzen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Darstellung von Proteinsequenzen . . . . . . . . . . . . . . . . . . . . . 9
2.4 Proteindomanen und -familien . . . . . . . . . . . . . . . . . . . . . . . 11
2.5 Datenbanken . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.5.1 SCOP (Structural Classification Of Proteins) . . . . . . . . . . . 14
2.5.2 Pfam (Protein Families) . . . . . . . . . . . . . . . . . . . . . . 16
3 Grundlagen des Sequenzalignments 17
3.1 Sequenzvergleich . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2 Ahnlichkeit und Distanz . . . . . . . . . . . . . . . . . . . . . . . . . . 20
iv
3.2.1 Hamming-Abstand und -Ahnlichkeit . . . . . . . . . . . . . . . 21
3.2.2 Levenshtein-Distanz . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.3 Komplexitat vs. Optimum . . . . . . . . . . . . . . . . . . . . . 24
3.3 Dynamisches Programmieren . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4 Paarweises Sequenzalignment . . . . . . . . . . . . . . . . . . . . . . . 28
3.4.1 Globales Alignment . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.4.2 Free-Shift Alignment . . . . . . . . . . . . . . . . . . . . . . . . 30
3.4.3 Lokales Alignment . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.5 Lucken und deren Bewertung . . . . . . . . . . . . . . . . . . . . . . . 32
3.6 Substitutionsmatrizen . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.6.1 PAM-Matrizen . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.6.2 BLOSUM-Matrizen . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.7 Scoring von Alignments . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.7.1 Empirische Festlegung eines Bezugssystems fur Scores . . . . . . 39
3.7.2 Verteilung der Scores . . . . . . . . . . . . . . . . . . . . . . . . 40
3.7.3 Standardisierung der Scores . . . . . . . . . . . . . . . . . . . . 42
4 Multiples Sequenzalignment 43
4.1 Anwendung multipler Sequenzalignments . . . . . . . . . . . . . . . . . 44
4.2 Komplexitat multipler Sequenzalignments . . . . . . . . . . . . . . . . 45
4.3 Iterative Generierung eines multiplen Alignments . . . . . . . . . . . . 47
4.4 Multiples Sequenzalignment mit ClustalW . . . . . . . . . . . . . . . . 48
4.5 Bewertung eines multiplen Sequenzalignments . . . . . . . . . . . . . . 50
5 Hidden Markov Modelle 51
5.1 Markov-Ketten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.2 Hidden Markov Modelle . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.3 Profil-HMMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.4 Verwendung von Profil-HMMs . . . . . . . . . . . . . . . . . . . . . . . 58
5.4.1 Training eines Profil-HMM . . . . . . . . . . . . . . . . . . . . . 58
5.4.2 Decodierung eines Profil-HMM . . . . . . . . . . . . . . . . . . 60
5.4.3 Alignment mit einem Profil-HMM . . . . . . . . . . . . . . . . . 60
v
6 Implementierung eines MSA mit einem Profil-HMM 62
6.1 Kurzbeschreibung der Implementierung . . . . . . . . . . . . . . . . . . 63
6.2 Generieren eines MSA mit einem Profil-HMM . . . . . . . . . . . . . . 63
6.2.1 Training des Profil-HMMs . . . . . . . . . . . . . . . . . . . . . 64
6.2.1.1 Berechnung der Emissionswahrscheinlichkeiten . . . . . 64
6.2.1.2 Berechnung der Ubergangswahrscheinlichkeiten . . . . 65
6.2.2 Single-Scoring und -Alignment mit dem Profil-HMM . . . . . . 66
6.2.3 Auswahl der zu alignierenden Sequenzen . . . . . . . . . . . . . 68
6.2.4 Progressives Alignment mit dem Profil-HMM . . . . . . . . . . 71
6.2.5 Kombination mehrerer Alignmentmethoden . . . . . . . . . . . 73
6.2.6 Progressive Expansion eines Profil-HMM . . . . . . . . . . . . . 75
7 Bewertung der Ergebnisse 79
7.1 Testdaten und -settings . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.2 Dichtefunktion der Scores . . . . . . . . . . . . . . . . . . . . . . . . . 82
7.3 Standardisierung der Dichtefunktion . . . . . . . . . . . . . . . . . . . 84
7.4 Drift der Scores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
7.5 Graphische Bewertung der Scores . . . . . . . . . . . . . . . . . . . . . 86
7.6 Numerische Bewertung der Scores . . . . . . . . . . . . . . . . . . . . . 90
7.7 Streuung der Scores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
7.8 Entropie der Emissionsmatrix . . . . . . . . . . . . . . . . . . . . . . . 94
7.9 Zusammenfassung der Evaluationsergebnisse . . . . . . . . . . . . . . . 99
8 Zusammenfassung und Ausblick 100
8.1 Zusammenfassung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
8.2 Ausblick . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
Literaturverzeichnis 103
Abkurzungsverzeichnis 108
Anhang 110
A Tabellen und Abbildungen 111
A.1 PAM250 Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
A.2 Score-Streudiagramme . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
vi
B Umgebung und Applikationen 115
B.1 Entwicklungs- und Testumgebung . . . . . . . . . . . . . . . . . . . . . 115
B.2 grepseq: Optionsbeschreibung . . . . . . . . . . . . . . . . . . . . . . . 116
B.3 amodseq: Optionsbeschreibung . . . . . . . . . . . . . . . . . . . . . . . 118
B.4 amodseq.ini: ClustalW Defaulteinstellungen . . . . . . . . . . . . . . . 119
C Daten- und Ergebnisdateien 120
C.1 Astral DB-Datei im FASTA Format . . . . . . . . . . . . . . . . . . . . 120
C.2 HMMOO Alignment Result Datei . . . . . . . . . . . . . . . . . . . . . 121
C.3 Z-Score Auswertung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
C.3.1 Z-Score Auswertung Zyklus 1 . . . . . . . . . . . . . . . . . . . 122
C.3.2 Z-Score Auswertung Zyklus 3 . . . . . . . . . . . . . . . . . . . 123
D Quelltexte 124
D.1 amodseq: Alignment voralignierter Sequenzen in das Profil-HMM . . . 124
E Datentrager 126
vii
Abbildungsverzeichnis
2.1 Darstellung der 3D-Struktur eines Proteins . . . . . . . . . . . . . . . . 6
2.2 Proteinsequenz eines Testdatensatzes im FASTA-Format . . . . . . . . 11
3.1 Einfaches Alignment von zwei Sequenzen . . . . . . . . . . . . . . . . . 17
3.2 Histogramm der Score-Verteilung randomisierter Sequenzen . . . . . . . 41
5.1 Erweiterter Zustandsgraph einer Markov-Kette . . . . . . . . . . . . . . 54
5.2 Schematische Darstellung eines Profil-HMM . . . . . . . . . . . . . . . 57
5.3 Einfaches MSA mit einem Profil-HMM . . . . . . . . . . . . . . . . . . 61
6.1 Verteilung der Reverse Corrected Scores . . . . . . . . . . . . . . . . . 69
6.2 Unvollstandiges MSA auf Basis eines Profil-HMM . . . . . . . . . . . . 73
6.3 MSA auf Basis eines Profil-HMM und ClustalW . . . . . . . . . . . . . 75
6.4 Symbolische Darstellung eines progressiven HMM-Expansion . . . . . . 76
7.1 Langenverteilung aller Sequenzen der Testdatenbank . . . . . . . . . . 81
7.2 Dichtefunktionen des RCS fur die Familie b.121.4.1 . . . . . . . . . . . 83
7.3 Standardisierte Dichtefunktionen des RCS . . . . . . . . . . . . . . . . 85
7.4 Dichtefunktionen des SCS und RCS fur e.3.1.1 nach dem Zyklus 1 . . . 87
7.5 Dichtefunktionen des SCS und RCS fur e.3.1.1 nach dem Zyklus 3 . . . 88
7.6 Dichtefunktionen des RCS fur a.138.1.1 nach Zyklus 1 . . . . . . . . . . 89
7.7 Dichtefunktionen des RCS fur a.138.1.1 nach Zyklus 3 . . . . . . . . . . 89
7.8 Darstellung der Scores in Streudiagrammen . . . . . . . . . . . . . . . . 92
7.9 Streuung der transformierten Scores nach einer PCA . . . . . . . . . . 93
7.10 Anderung der transformierten Scores eines Profil-MSA . . . . . . . . . 97
viii
Tabellenverzeichnis
2.1 Standardisiertes Alphabet zur Codierung der Proteinsequenzen . . . . . 10
3.1 Matrix zur Berechnung des Levenshtein-Abstands . . . . . . . . . . . . 26
3.2 Auszug aus einer PAM250 Matrix . . . . . . . . . . . . . . . . . . . . . 35
3.3 BLOSUM-62 Substitutionsmatrix . . . . . . . . . . . . . . . . . . . . . 38
7.1 Multiple Alignments zur Beschreibung von Profil-HMMs . . . . . . . . 80
7.2 Globale Mittelwerte der RCS und SCS nach 4 Zyklen . . . . . . . . . . 91
7.3 Vergleich der Dichte- mit der Entropieanderung . . . . . . . . . . . . . 98
A.1 Vollstandige PAM250 Matrix . . . . . . . . . . . . . . . . . . . . . . . 112
A.2 Scores der Familie a.1.1.2 nach dem ersten Zyklus . . . . . . . . . . . . 113
A.3 Scores der Familie a.1.1.2 nach dem dritten Zyklus . . . . . . . . . . . 114
ix
1
Einfuhrung
Komplexe Biomolekule spielen im Bauplan des Lebens eine entscheidende Rolle. Sol-
che Makromolekule setzen sich dabei typischerweise aus einer Sequenz von Bausteinen
zusammen. Typische Vertreter aus der Biologie sind die in allen lebenden Organis-
men vorkommenden und deshalb als naturliche Nukleinsauren bezeichnete DNA1 und
RNA2 oder die darauf basierenden Eiweißstoffe (Proteine) aus einer Abfolge von Ami-
nosauren. In der Molekularbiologie werden die Teile dieser Biomakromolekule als Folge
von Symbolen beschrieben, wobei eine Sequenz dieser Symbole eine lineare Abfolge von
Basen oder Aminosauren definiert. In der Bioinformatik werden diese Sequenzen zur
Verarbeitung in Form einfacher Zeichenketten (Strings) kodiert, welche eine wesent-
liche Datengrundlage aller bioinformatischen Sequenzanalysen und -interpretationen
darstellen.
Der Aufbau eines Makromolekuls und die damit einhergehende Abfolge bestimmter
Symbole in den Sequenzen werden von Biologen oft auch fur bestimmte Eigenschaften
und Funktionen verantwortlich gemacht. Heute weiß man, dass eine ahnliche Abfolge
bestimmter Symbole in einer Sequenz oder in Sequenzteilen auf ahnliche Eigenschaf-
ten und Funktionen des Makromolekuls hinweisen kann [32]. Eines der wesentlichen
Anliegen der Bioinformatik ist deshalb die Klarung der Frage, ob eine bestimmte Se-
quenz einer oder mehreren anderen Sequenzen ahnlich ist. Findet man durch den Ver-
gleich mehrerer verwandter Proteinsequenzen ahnliche Organisationseinheiten, deren
1Desoxyribonukleinsaure2Ribonukleinsaure
1
1. Einfuhrung 2
Kombination die Eigenschaften eines Proteins maßgeblich bestimmt, so konnen ver-
wandte Sequenzen dadurch bestimmten Familien und Funktionen zugeordnet werden.
Zu beantworten ist dabei immer die Frage, was Ahnlichkeit in diesem Zusammenhang
uberhaupt heißen kann.
1.1 Problemstellung und Motivation
Funktionelle Ahnlichkeiten verschiedener Sequenzen basieren meist auf evolutionaren
Verwandtschaften (Homologien) der Sequenzen. Die evolutionaren Entwicklungen ha-
ben im Laufe der Zeit eventuell eine Sequenz verandert, dies hat aber nicht zwangsweise
eine wesentliche Anderung der grundlegenden Eigenschaften und Funktionen der mu-
tierten Proteine zur Folge. Fur die Suche nach Regelmaßigkeiten und Ahnlichkeiten
werden vorzugsweise statistische Methoden eingesetzt, welche das Auftreten von Sym-
bolen in einer oder mehreren Sequenzen und -positionen bewerten. Zusatzlich wird
festgestellt, auf welche Sequenzteile andere Teile folgen, oder ob einzelne Sequenz-
bausteine durch andere ersetzt, entfernt oder eingefugt wurden. In der Bioinformatik
versucht man nun, die Ahnlichkeit oder den”Abstand“ von Sequenzen numerisch aus-
zudrucken. Kleinere Abstande weisen dabei auf großere Ahnlichkeiten und eine evolu-
tionare Nahe und großere Abstande auf großere Unahnlichkeiten hin. Dadurch werden
Sequenzen untereinander vergleichbar und Ahnlichkeiten erkennbar, selbst wenn die
einzelnen Bausteine fur sich und vorweg betrachtet nicht ahnlich erscheinen wurden.
Ein Alignment von Sequenzen zielt nun darauf ab, die in den Sequenzen vorkommenden
Muster zu identifizieren und deren Position innerhalb der Sequenzen zu bestimmen.
Beim Sequenzalignment wird versucht, zwei oder mehrere Bausteinketten aufgrund
ihrer Ahnlichkeit zueinander auszurichten. Großere Ahnlichkeiten innerhalb von Se-
quenzen weisen auf großere Gemeinsamkeiten der DNA-Muster hin, wodurch sich auch
großere Bereiche zueinander gut alignieren lassen. Durch das Alignment sollen gleiche
oder ahnliche Teile der Strings bzw. Regionen gleicher Funktionalitat untereinander no-
tiert werden konnen. Die Reihenfolge der Symbole muss dabei erhalten bleiben, wobei
- den Mutationsformen der Natur folgend - auch Leerstellen (Gaps) in die Sequenzen
1. Einfuhrung 3
eingefugt werden konnen, um Einfugungen in anderen Sequenzen auszugleichen und
eine spaltenweise Alignierung zusammenhangender Sequenzabschnitte zu ermoglichen.
Das Alignment von Sequenzen gilt im Allgemeinen als nicht-triviale Problemstellung.
Der paarweise Sequenzvergleich, also der Vergleich von nur zwei Sequenzen, kann schon
einen betrachtlichen Aufwand generieren, weist aber eine noch vergleichsweise niedrige
algorithmische Komplexitat auf. Die Alignierung ganzer Sequenzgruppen fuhrt zu deut-
lich besseren Ergebnissen und einer hoheren Sensitivitat der Algorithmen. Der Aufwand
dafur ist jedoch schon bedeutend großer, denn hinter einem effizienten und zielsicheren
Alignmentverfahren stecken meist komplexe numerische, algorithmische und mathema-
tische Methoden und Uberlegungen, welche die Forschung in der Bioinformatik pragen.
Das Ziel aller Uberlegungen ist einerseits eine performante Alignierungsmethode zu ent-
wickeln, andererseits eine maximale Ubereinstimmung der alignierten Sequenzen unter
der Pramisse eines biologisch bzw. evolutionar sinnvollen Alignments zu erreichen.
1.2 Uberblick
Das folgende Kapitel 2 fuhrt in ausgewahlte Grundlagen der Bioinformatik ein
und erklart die zum Verstandnis notwendigen Begriffe sowie den Aufbau von Proteinen.
Nachdem Proteinsequenzen in der Bioinformatik ublicherweise in Form von Zeichenket-
ten verarbeitet werden, nimmt die Codierung und Speicherung von Proteinsequenzen
einen wesentlichen Teil dieses Kapitels ein.
In Kapitel 3 werden die Grundlagen des Sequenzalignments vorgestellt und die Be-
griffe Ahnlichkeit und Distanz diskutiert. Einfache Beispiele fuhren in die grundlegen-
den Algorithmen des paarweisen Sequenzalignments ein und zeigen einfache Scoring-
Methoden zur objektiven und quantitativen Bewertung von Ahnlichkeit und Distanz.
Das Kapitel 4 widmet sich den multiplen Sequenzaligments (MSA) und der Kom-
plexitat multipler Alignments. Am Beispiel ClustalW wird eine iterative Methode vor-
gestellt, in der Sequenzen schrittweise zu einem MSA aligniert werden.
1. Einfuhrung 4
Im Kapitel 5 werden die mathematischen Grundlagen eines Hidden Markov Modells
(HMM) geklart. Uber Markov-Ketten werden die Begriffe Emissions- und Translations-
wahrscheinlichkeit nahergebracht, bevor die schrittweise Anwendung eines Profil-HMM
vom Training uber die Decodierung bis zum Alignment behandelt wird.
Das Kapitel 6 stellt die Implementierung eines MSA auf Basis eines Profil-HMM
vor. Mit Bezug auf die theoretischen Grundlagen werden die Einzelschritte zur Er-
stellung eines multiplen Alignments erlautert und die entwickelten Algorithmen und
Applikationen beschrieben.
Kapitel 7 widmet sich der Bewertung der generierten Alignments. Der Schwerpunkt
liegt dabei nicht in der Prufung der biologischen Qualitat der Applikationen, sondern
in erster Linie in der Interpretation der Scores zur Bestimmung der Qualitat der Align-
ments und in der Beeinflussung des Modells durch die Anderung des Profil-HMM in
einem iterativen Prozess.
Das letzte Kapitel 8 schließt die Arbeit mit einer Zusammenfassung ab, die einen
kurzen Uberblick uber die entwickelten Methoden und Applikationen, die Bewertung
der Ergebnisse und die gewonnenen Erkenntnisse gibt. Im Ausblick werden mogliche
Erweiterungen und Verbesserungen der im Rahmen dieser Arbeit entwickelten Metho-
den und Applikationen vorgeschlagen.
2
Allgemeine Grundlagen
Bioinformatik ist die Entwicklung und das Betreiben von Datenbanken,
Software und mathematischen Werkzeugen zur Analyse, Organisation und
Interpretation biologischer Daten. [21, S. 3]
Das Ziel dieser Arbeit ist die Entwicklung eines Verfahrens zur Suche von strukturel-
len Ahnlichkeiten zweier oder mehrere Proteine. Die Ahnlichkeiten werden auf Basis
von Sequenzen ermittelt und stellen ein Abbild der funktionellen Verwandtschaft und
genetischen Entwicklung dar. Auf Grundlage der Ahnlichkeiten werden die Sequenzen
spaltenweise so untereinander ausgerichtet, dass jene Bereiche sichtbar werden, in der
- vorerst nur hypothetisch - sowohl strukturelle als auch funktionelle Ahnlichkeiten zu
erwarten sind [32]. Bevor auf die Details eines Alignments naher eingegangen werden
kann, mussen einige Grundlagen erlautert und Begriffe definiert werden.
Im aktuellen Kapitel werden - Hutt und Dehnert’s Definition von Bioinformatik fol-
gend - allgemeine bioinformatische Grundlagen eingefuhrt, soweit sie zum Verstandnis
und der Darlegung der Problemstellung dienlich sind. Auf molekularbiologische und
biochemische Details wird dabei mit dem Fokus auf das Wesentliche bewusst verzich-
tet. Sie sind bei Bedarf in etwaigen fachspezifischen Quellen nachzulesen. Geklart wird
im Folgenden vielmehr, wie bioinformatische Sequenzen dargestellt und gespeichert
werden, welche Informationen aus ihnen gewonnen werden konnen und was man unter
einem Sequenzvergleich versteht.
5
2. Allgemeine Grundlagen 6
2.1 Proteinstruktur
Im Wesentlichen werden die Eigenschaften eines Makromolekuls durch dessen dreidi-
mensionale Struktur (Raumstruktur) bestimmt. Die Raumstruktur ist wiederum die
Folge einer Sequenz aus unterschiedlichen Grundbausteinen und deren Anordnung der
Atome im Raum.
Bei Proteinen beschreibt die Abfolge der unterschiedlichen Aminosauren die sogenann-
te Primarstruktur des Makromolekuls, wobei bestimmte Strukturteile von Biologen fur
bestimmte Funktionen verantwortlich gemacht werden. Aus einer bestimmten Anord-
nung der Aminosauren innerhalb einer Proteinsequenz entstehen durch den Prozess
der Faltung in der Regel eindeutig bestimmte dreidimensionale Substrukturen. Diese
Regionen bestimmen die elementaren dreidimensionalen Grundformen eines Proteins.
Abbildung 2.1: Darstellung der 3D-Struktur eines Proteins.
Typische strukturelle Grundelemente und Grundformen bei Proteinen sind die soge-
nannten α-Helices mit ihren spiralformigen Strukturen und β-Faltblatter, welche die
Faltung der dreidimensionalen Struktur im Raum ermoglichen. Man bezeichnet die-
se lokalen Strukturelemente als die Sekundarstruktur eines Proteins (siehe Abbildung
2.1). Zusammen mit den als Verbindungselemente fungierenden Zwischenregionen bil-
den sie die Tertiarstruktur eines Proteins, welche die raumliche Anordnung der Atome
beschreibt. Die Verschrankung mehrerer Tertiarstrukturen zu einem Proteinkomplex
bildet die Quaternarstruktur oder Quartarstruktur [32, 21].
2. Allgemeine Grundlagen 7
Versuchen Molekularbiologen nun den Aufbau und die Funktion eines Gens oder Pro-
teins zu entschlusseln, so weichen sie vorzugsweise auf die Untersuchung der Sequen-
zen aus. Die Ermittlung einer Sequenz als Ausgangsbasis weiterer Untersuchungen ist
vergleichsweise kostengunstig, stellt man sie der Vermessung der Raumstruktur mit
aufwendigen und deutlich teureren Methoden gegenuber. Zudem weisen homologe Se-
quenzen, also ahnliche bzw. verwandte Gene mit gleichen Vorfahren, oft auch ahnliche
3D-Strukturen und eine ahnliche Anordnung von Symbolen auf. Merkl und Waack
beschreiben dieses Zusammenhang kurz als den Generellen Grundsatz: Sequenz deter-
miniert Struktur. [32, S. 153]
2.2 Sequenzen
In der Molekularbiologie werden DNA-Sequenzen oder Proteine als Folge von Sym-
bolen beschrieben. Eine Sequenz solcher Symbole definiert dabei eine lineare Abfolge
kleinerer Molekulbausteine, bestehend aus Basen oder Aminosauren. Im Laufe der Evo-
lution werden diese Sequenzen variiert, indem Teile der Molekulstrukturen und damit
auch der Sequenzen dupliziert oder modifiziert werden. Oft werden in sogenannten
Punktmutationen einzelne Symbole, oft ganze Sequenzfragmente gegen andere ausge-
tauscht oder innerhalb der Sequenz verschoben. Gleiche oder ahnliche Strukturen oder
evolutionar”erfolgreiche“ Bausteine werden oft auch wiederverwendet und tauchen in
unterschiedlichsten biologischen Sequenzen wiederholt auf [32].
Mittlerweile gehoren Sequenzen zum wichtigsten Datenmaterial der Bioinformatik. Ein
Grund dafur sind die in den Siebzigerjahren des letzten Jahrhunderts von Gilbert [16]
und Sanger [40] entwickelten Techniken1, die ein automatisiertes, schnelles und vor al-
lem aber kostengunstiges Sequenzieren von DNA erst moglich gemacht haben [8]. Seit
etwa Mitte der 1980er Jahre werden automatische Sequenzierer kommerziell eingesetzt.
1Der zu dieser Zeit auf der Harvard University arbeitende Biologe Walter Gilbert und der inCambridge forschende Molekularbiologe Frederick Sanger haben
”for their contributions concerning
the determination of base sequences in nucleic acids“ [36] im Jahr 1980 einen Nobelpreis fur Chemieerhalten.
2. Allgemeine Grundlagen 8
Heute werden mit dem sogenannten Shotgun-Sequencing meistens zufallige Teilsegmen-
te einer langeren DNA-Sequenz sequenziert und dann beim sogenannten Sequence As-
sembling die uberlappenden Einzelteile mit speziellen Algorithmen zusammengefugt
[21, 32]. Um Sequenzen softwaretechnisch verarbeiten und austauschen zu konnen, be-
darf es einer Abstraktion der Sequenzen und speziellen Codierung der Grundbausteine.
Eine DNA oder RNA als Sequenz dargestellt besteht aus einer Abfolge von vier Ba-
sen. Diese hintereinander gereiht bilden den genetischen Bauplan bzw. im Falle eines
Genoms die komplette Erbinformation eines Lebewesens. Den Bioinformatikern dient
eine Abfolge von Buchstaben der Abstraktion von Sequenzen. DNA-Sequenzen werden
beispielsweise aus den vier Zeichen A, C, G und T gebildet, welche fur vier Basen2
stehen. Die Menge ΣD der Zeichen zur Beschreibung der DNA und die Menge ΣR der
Zeichen zur Beschreibung der RNA lautet:
ΣD = {A,G,C,T} ΣR = {A,G,C,U}
Unter Transkription bezeichnet man allgemein die Ubertragung eines Textes von einem
System in ein anderes System. In Analogie dazu wird der Ubersetzungsprozess von einer
DNA-Sequenz in eine RNA-Sequenz von Biologen ebenso Transkription genannt [32].
Aus der Sicht der Symbole entspricht eine Transkription einer DNA- in eine RNA-
Sequenz formal einem Alphabetwechsel von ΣD nach ΣR [21]. Vereinfacht dargestellt
entsteht aus einer DNA-Sequenz durch den Austausch der Base T (Thymin) gegen eine
Base U (Uracil) eine RNA-Sequenz.
ΣD = {A,G,C, T} A G T C T C G T T A C T T C T T C A
T −→ U Transkription
ΣR = {A,G,C, U} A G U C U C G U U A C U U C U U C A
Das hierbei entstehende RNA-Molekul wird Transkript genannt [32] und tragt zur
Ubersetzung der genetischen Informationen in Proteine bei.
Bei Proteinen bilden Aminosauren die Grundbausteine der Proteinsequenzen. Die Dar-
stellung von Proteinsequenzen erfolgt mittels eines standardisierten Alphabets mit dem
sogenannten One-Letter-Code.
2Die Buchstaben zur Darstellung einer DNA- oder RNA-Sequenz stehen dabei jeweils fur die Namenvon Basen. Dabei steht fur Adenin ein ’A’, fur Guanin ein ’G’, fur Thymin ein ’T’, fur Cytosin ein’C’ und fur Uracil ein ’U’.
2. Allgemeine Grundlagen 9
2.3 Darstellung von Proteinsequenzen
Aus der Sicht der Bioinformatiker basieren auch Untersuchungen von Proteinen oft auf
den Analysen und Interpretationen einfacher Zeichenketten. Sie stellen die abstrahier-
ten Sequenzen realer Aminosaureketten dar.
Auch Proteine sind Makromolekule. Diese werden aus naturlich vorkommenden Ami-
nosauren gebildet, welche uber den genetischen Code durch jeweils drei Basen einer
DNA bestimmt werden. Eine derartige Dreierkombination wird als Triplett oder Co-
don bezeichnet. Aus den vier moglichen Zeichen wurden sich 43 = 64 Dreierkombi-
nationen bilden lassen, aufgrund biologischer Gegebenheiten, auf die hier nicht naher
eingegangen wird, definieren die 64 moglichen Tripletts jedoch nur 20 unterschiedliche
Aminosauren [32].
Besonders beim Sequenzalignment erweisen sich Tripletts als nicht sonderlich gunstig,
da die ublicherweise zur Anwendung kommenden Verfahren die Sequenzen zeichen-
weise abarbeiten. Alignierte Sequenzen werden an bestimmten Stellen oft mit Lucken
(Gaps) aufgefullt, um eine bessere Alignierung zu erreichen. Damit wurden gerade
beim Alignment von Proteinsequenzen die Tripletts der einzelnen Aminosauren auf-
gebrochen. Weist man jeder dieser Aminosauren aber nur ein bestimmtes Zeichen der
Menge ΣA zu, so bleiben die Dreierkombinationen der jeweiligen Tripletts auch bei ei-
nem zeichenweisen Alignment bestehen. Die Menge der Zeichen zur Beschreibung von
Aminosauresequenzen ist standardisiert als [22]:
ΣA = {A,C,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y}
Dieser Wechsel des Alphabets wird als Translation bezeichnet. Jeder Buchstabencode
der Menge ΣA steht fur eine bestimmte Aminosaure.
ΣD = {A,G,C, T} A G T︸ ︷︷ ︸ C T C︸ ︷︷ ︸ G T T︸ ︷︷ ︸ A C T︸ ︷︷ ︸ T C T︸ ︷︷ ︸ T C A︸ ︷︷ ︸Triplett−→Aminosaure Translation
ΣA = {A,C . . . ,W, Y } S L V T S S
Eine vollstandige Liste mit den Buchstabencodes, den jeweiligen Aminosauren und
deren Codone zeigt die Tabelle 2.1.
2. Allgemeine Grundlagen 10
Buchstabencode Aminosaure Codone
A Alanin (Ala) GCA, GCC, GCG, GCT
C Cystein (Cys) TGC, TGT
D Aspartat (Asp) GAC, GAT
E Glutamat (Glu) GAA, GAG
F Phenylalanin (Phe) TTC, TTT
G Glycon (Gly) GGA, GGC, GGG, GGT
H Histidin (His) CAC, CAT
I Isoleucin (Ile) ATA, ATC, ATT
K Lysin (Lys) AAA, AAG
L Leucin (Leu) CTA, CTC, CTG, CTT, TTA, TTG
M Methionin (Met) ATG
N Asparagin (Asn) AAC, AAT
P Prolin (Pro) CCA, CCC, CCG, CCT
Q Glutamin (Gln) CAA, CAG
R Arginin (Arg) AGA, AGG, CGA, CGC, CGG, CGT
S Serin (Ser) AGC, AGT, TCA, TCC, TCG, TCT
T Threonin (Thr) ACA, ACC, ACG, ACT
V Valin (Val) GTA, GTC, GTG, GTT
W Tryptophan (Trp) TGG
Y Tyrosin (Tyr) TAC, TAT
Tabelle 2.1: Standardisiertes Alphabet zur Codierung der Proteinsequenzen [22].
In Tabelle 2.1 werden die Dreierkombinationen (Codone) und die jeweiligen Ami-
nosauren, die daraus gebildet werden, gelistet. Betrachtet man die Codonspalte der
Aminosauren so fallt auf, dass bestimmte Aminosauren aus unterschiedlichen Codone
bestehen konnen. Ebenso auffallend ist, dass eine Mutation der dritten Base im Codon
die Art der Aminosaure nicht zwingend andert. Die großten Anderungen der Eigen-
schaften erfahrt eine Aminosaure, wenn die mittlere der drei Basen mutiert. Die in
der Tabelle nicht aufgelisteten Tripletts sind keine Aminosauren, sondern ubernehmen
andere Funktionen innerhalb einer DNA-Sequenz. Als Beispiele konnen spezielle Co-
done angefuhrt werden, die den Beginn (Startcodon) oder das Ende (Stoppcodon) einer
Proteinsequenz markieren [32].
2. Allgemeine Grundlagen 11
Mit diesem One-Letter-Code kann ein Protein mit einer auf ein Drittel der ursprung-
lichen Lange verkurzten Zeichenkette (String) eindeutig beschrieben werden. Die Ab-
bildung 2.2 zeigt die Sequenz des Proteins d1ngka in der 1-Zeichen-Codierung:
>d1ngka_ a.1.1.1 (A:) Protozoan/bacterial hemoglobin
ksfydavggaktfdaivsrfyaqvaedevlrrvypeddlagaeerlrmfleqywggprty
seqrghprlrmrhapfrislierdawlrcmhtavasidsetlddehrrelldylemaahs
Abbildung 2.2: Proteinsequenz eines Testdatensatzes im FASTA-Format. In der erstenZeile markiert das Zeichen ’>’ den Beginn einer Sequenz mit dem Bezeichner, ab derzweiten Zeile beginnt die Beschreibung der Sequenz im One-Letter-Code.
Wenn im Folgenden von Sequenzen geschrieben wird, dann sind, wenn nicht aus-
drucklich anders vermerkt, immer Proteinsequenzen gemeint, auf deren Analyse der
Schwerpunkt dieser Arbeit liegt.
So zufallig und chaotisch Proteine und ihre Sequenzen auf den ersten Blick erscheinen
mogen, so folgen sie in ihrer Struktur und ihrer Symbolabfolge doch einem System.
Im direkten Vergleich weisen Sequenzen in Abhangigkeit ihrer Entstehungsgeschichte
Bereiche mit strukturellen Regelmaßigkeiten und Gemeinsamkeiten auf, die eine Ein-
teilung in Proteindomanen und Proteinfamilien nahelegen.
2.4 Proteindomanen und -familien
Vergleicht man die 3D-Struktur zweier Proteine, so zeigt sich, dass bestimmte Bereiche
eine großere Ahnlichkeit aufweisen, in anderen Bereichen wiederum nur geringe oder
keinerlei Ahnlichkeiten erkennbar sind. Diese Ahnlichkeitsschwankungen (im Bereich
von signifikant hohen Ahnlichkeitsgraden bis hin zu volliger Verschiedenheit) werden
auch beim Sequenzvergleich sichtbar. Dabei variiert die Ahnlichkeit stark uber die
Lange der Sequenzen. Der Grund dieser partiellen Schwankungen liegt im modularen
Aufbau eines Proteins, welches sich aus großeren Einheiten zusammensetzt. Die klein-
ste Einheit einer definierten und geometrisch relativ unabhangigen gefalteten Struk-
tur bestehend aus α-Helices und β-Faltblatter nennt man strukturelle Domane. Sind
bestimmte Sequenzabschnitte fur die Eigenschaften eines Proteins verantwortlich, so
wird sie als funktionelle Domane bezeichnet. Strukturelle und funktionelle Domanen
2. Allgemeine Grundlagen 12
mussen sich nicht zwangsweise decken, auch wenn die Erfahrung zeigt, dass dies oft
der Fall ist. Die meist aus etwa 50 bis 150 Aminosauren bestehenden Domanen bilden
die Grundbausteine großerer Proteine. Im Laufe der evolutionaren Entwicklung werden
Proteindomanen haufig wiederverwendet und zu neuen Proteinen kombiniert [32].
Unter Berucksichtigung der Kombinationsmoglichkeiten der 20 unterschiedlichen Ami-
nosauren liegt der Gedanke nahe, dass es in der Natur eine nahezu unerschopfliche
Proteinvielfalt geben musste. In der Praxis zeigt sich aber, dass sich bestimmte Kom-
binationen im Laufe der Evolution mehr, andere nicht durchsetzen konnten. Die An-
zahl der Grundbausteine der Natur ist stark eingeschrankt. Ein Großteil der bekannten
Proteine besteht aus einer Domane (Eindomanenproteine) und schatzungsweise 80 %
aller Proteine werden einem von 400 Faltungstypen zugeordnet [32]. Klassifiziert man
Proteine nach deren gemeinsamen Vorfahren, so konnen sie bestimmten Proteinfa-
milien zugerechnet werden. Proteine gleicher Familien haben in der Regel ahnliche
3D-Strukturen und Funktionen und weisen eine signifikante Sequenzahnlichkeit auf.
Trotz aller Einschrankungen der Natur ist die Anzahl verschiedener Sequenzen und
Domanenkombinationen betrachtlich und die Organisation der Datenbestande aufwen-
dig. Sequenzen werden deshalb nach verschiedenen Kriterien, wie beispielsweise Prote-
indomanen und Proteinfamilien katalogisiert, und zum großten Teil, wenn keine patent-
rechtlichen oder wirtschaftlichen Grunde dagegensprechen, in offentlich zuganglichen
Datenbanken abgelegt.
2.5 Datenbanken
Im Laufe der Jahrzehnte wurde eine Menge biologischer und biochemischer Daten ge-
neriert, die in verschiedensten Datenbanken abgelegt wurden. Diese decken mit ihrem
Angebot die speziellen Bedurfnisse unterschiedlichster Forschungsbereiche ab, wobei
den Sequenz-Datenbanken eine besondere Bedeutung zukommt. Gespeichert werden in
diesen Datenbanken nicht nur die Daten sequenzierter Gene, sondern auch zusatzliche
(Meta-)daten, welche je nach Anwendungsfeld fur die Auswertung, die Suche in den Da-
tenbanken und den Vergleich der Datensatze untereinander verwendet werden konnen.
2. Allgemeine Grundlagen 13
Mittlerweile nehmen die Datenbestande einen nicht unbetrachtlichen Umfang an, dem-
entsprechend viel Aufwand wird datenbankseitig investiert und dementsprechend viel-
seitig ist das Angebot. Die Große einiger Datenbanken bewegt sich mittlerweile in Rich-
tung Petabyte3 [3]. Ob ihrer Anzahl und Vielseitigkeit kann man schnell den Uberblick
verlieren. So stellen beispielsweise Cochrane und Galperin im Dezember 2009 in [9]
fest, dass mit Ende des Jahres insgesamt 1230 Datenbanken in der Database Issue and
Database Collection registriert4 sind.
Fur statistische Auswertungen ist die Menge der zur Verfugung stehenden Datensatze
wesentlich. Eine Vernetzung der Datenbestande untereinander sichert die Qualitat der
Daten und lasst komplexe Analysen zu, jedoch kommt neben der Datenmenge damit
der einheitlichen Auszeichnung der Datenbestande wie auch der Qualitat der Daten
eine besondere Bedeutung zu [32].
Aufgrund der rasant zunehmenden Datenbankgroßen ist es nicht verwunderlich, dass
immer wieder neue Methoden zur effizienten Speicherung und Filterung der Daten ent-
wickelt werden. Herkommliche Abfragesprachen wie SQL (Structured Query Language)
finden selten Anwendung, wenn aus der Menge aller Sequenzen nur bestimmte Sequen-
zen gefunden werden sollen oder Sequenzen untereinander verglichen werden mussen.
Die Filter verfolgen vielmehr das Prinzip der Suche nach Ahnlichkeiten, welche in ir-
gendeiner Form quantifiziert werden mussen (siehe auch Kapitel 3.1 und Kapitel 3.2ff).
FASTA5 und BLAST6, die beiden bekanntesten Sequenzfilter, vergleichen bei einer Ab-
frage beispielsweise samtliche Eintrage der Datenbank mit der Ausgangssequenz (query
sequence) und errechnen die Ahnlichkeiten mithilfe von Scoring-Matrizen, auch Substi-
tutionsmatrizen (siehe Kapitel 3.6) genannt, um die Menge der Zielsequenzen (target
sequences) zu ermitteln [21, 31, 33]. Beide Algorithmen haben gemeinsam, dass sie
in einem Vorauswahlverfahren mittels Approximation aus allen Sequenzen eine stark
3Ein Petabyte sind 1000 Terabyte oder anders ausgedruckt etwa 10 hoch 15 Byte.4Alljahrlich wird im ersten Heft der Zeitschrift Nucleic Acids Research eines jeden Jahrgangs in der
”Database Issue and Database Collection“ ein Uberblick uber die Entwicklung molekularbiologischer
Datenbanken gegeben [9]. Eine aktuelle Liste der gesammelten Datenbanken findet man unter OnlineDatabase Collection, http://www.oxfordjournals.org/nar/database/a/ (Stand: 7. April 2010).
5FASTA (FAST-ALL) ist eine Suchmethode zur schnellen Filterung von Sequenzen in Datenbanken:http://www.ebi.ac.uk/Tools/fasta/ (Stand: 3. August 2010)
6BLAST (Basic Local Alignment Search Tool): http://www.ebi.ac.uk/Tools/blast/ (Stand: 3.August 2010)
2. Allgemeine Grundlagen 14
verkleinerte Menge der potentiell interessanten Sequenzen ermittelt, welche dann ei-
ner genaueren Analyse unterzogen werden. Ohne diese Vorauswahl waren On-the-Fly-
Abfragen in großen Datenbanken kaum mehr moglich [33].
Bei den Datenbanken unterscheidet man den abgelegten Daten zufolge zwischen prima-
ren und sekundaren Datenbanken. In den primaren Datenbanken werden die experi-
mentell ermittelten Rohdaten abgelegt, wohingegen in den sekundaren Datenbanken
das von den Sequenzen abgeleitete Wissen samt allen Querverweisen gespeichert wird.
Eine der erfolgreichsten Protein-Datenbanken dieser Art ist SWISS-PROT7, welche
in ihren Datensatzen wiederum auf mehr als 100 andere Datensammlungen referen-
ziert [46, 32]. Auf Basis der Spezialisierung auf bestimmte Daten bzw. der Darstellung
der Sequenzen werden sie auch in DNA-Sequenz-, RNA-Sequenz-, Proteinsequenz und
Proteinstrukturdatenbanken unterteilt.
Aufgrund der Fulle der Datenbanken kann hier nur noch einmal auf die Database Issue
and Database Collection8 von Nucleic Acids Research verwiesen werden. Exemplarisch
werden zwei Datenbanken herausgegriffen und kurz eingefuhrt, die im Rahmen dieser
Arbeit zur Anwendung kommen, SCOP und Pfam.
2.5.1 SCOP (Structural Classification Of Proteins)
Die Proteindatenbank SCOP9 (Structural Classification Of Proteins) organisiert ih-
re Daten anhand einer hierarchischen Klassifikation der Proteine basierend auf Se-
quenzahnlichkeiten und der 3D-Struktur eines Proteins. Ziel von SCOP ist die Abbil-
dung struktureller und evolutionarer Verwandtschaftsbeziehungen von Proteinen [2].
Im Laufe der evolutionaren Entwicklungen mutieren Proteine, wobei sich deren Se-
quenzen und Funktionen deutlich verandern konnen. Die Struktur der Proteine ist
oft starker konserviert als die Sequenzen und Funktionen es sind, und erlaubt einen
Ruckschluss auf gemeinsame evolutionare Entwicklungen, selbst wenn der Sequenzver-
gleich keine Verwandtschaftsbeziehungen mehr erkennen lasst. Murzin et al. [2] nutzten
7SWISS-PROT Protein Knowledgebase: http://expasy.org/sprot/ (Stand: 14. August 2010)8Database Issue and Database Collection: http://www.oxfordjournals.org/nar/database/a/
(Stand: 7. April 2010)9SCOP (Structural Classification of Proteins) Database (1.75 release, June 2009): http://scop.
mrc-lmb.cam.ac.uk/scop/ (Stand: 1. Juli 2010)
2. Allgemeine Grundlagen 15
diese Kenntnisse und klassifizierten 1995 Proteine auf Basis von Proteindomanen (sie-
he Abschnitt 2.4). Die Autoren fuhren in ein hierarchisches Klassifikationsschema ein,
welches vier Ebenen beschreibt. Diese Ebenen sind:
• Proteinfamilie (Family): In den Proteinfamilien werden jene Proteine zusammen-
gefasst, deren Sequenzen zu mindestens 30% ident sind oder solche, deren Se-
quenzahnlichkeiten geringer sind, jedoch deren Funktionen und Strukturen große
Ahnlichkeiten aufweisen.
• Superfamilie (Superfamily): Proteine, die geringe Sequenzahnlichkeiten aufwei-
sen, deren evolutionare Verwandtschaft aufgrund ahnlicher Funktionen und Struk-
turen jedoch wahrscheinlich ist, werden zu Superfamilien zusammengefasst.
• Faltungstyp (Common Fold): Familien und Superfamilien werden dann unter
einem Faltungstyp zusammengefasst, wenn die Anordnung der wesentlichen Se-
kundarstrukturelemente große Ahnlichkeiten und vergleichbare topologische Ver-
bindungselemente aufweisen.
• Klassen (Class): Faltungstypen werden in funf Klassen aufgeteilt, die sich aus
der Zusammensetzung der Sekundarstrukturelemente definieren. Hierbei unter-
scheidet man die Zusammensetzung der Proteinstrukturen nach: (a) einer all-
alpha (hauptsachlich aus α-Helices bestehend); (b) all-beta (hauptsachlich aus
β-Faltblatter bestehend); (c) alpha and beta (aus α-Helices und β-Faltblatter
bestehend); (d) alpha plus beta (aus mehrheitlich isolierten α-Helices und β-
Faltblattern bestehend); (e) multi-domain (fur Domanen mit unterschiedlichen
Faltungen und wenigen bekannten Homologien).
Fur die Beschreibung der SCOP-Hierarchie wurde eine spezielle Notation definiert.
Dabei wird die SCOP Hierarchie in einem Set of Concise Classificaton Strings (SCCS)
wie folgt beschrieben [4]:
class.fold.superfamily.family
Der SCCS ist eine kompakte Notation der SCOP Domanenklassifikation zur hierarchi-
schen Beschreibung der Klasse, des Faltungstyps, der Superfamilie und der Familie. Die
in Abbildung 2.2 dargestellte Sequenz des Proteins Protozoan/bacterial hemoglobin
weist einen SCCS von a.1.1.1 aus. Das Zeichen ’a’ beschreibt die Klasse (all alpha
2. Allgemeine Grundlagen 16
proteins), ’1’ den Faltungstyp (globin-like), ’1’ die Superfamilie (globin-like) und ’1’
die Familie (truncated hemoglobin), denen diese Sequenz zugeordnet wird.
Weitere auf SCOP basierte Datenbanken verwenden ebenfalls diese Notation. ASTRAL10
und die im Folgenden beschriebene Pfam sind nur zwei davon.
2.5.2 Pfam (Protein Families)
Pfam11 (Protein Families), eine 1997 von Sonnhammer et al. [44] eingefuhrte Proteinda-
tenbank, enthalt in der Version 24.0 eine Sammlung von 11912 Proteinfamilien, die als
multiple Sequenzalignments oder in Form von Hidden Markov Modellen reprasentiert
werden. Pfam besteht aus zwei Komponenten: Pfam-A liefert die Sequenzen von Fami-
lien, die manuell zugewiesen und aligniert wurden und jeweils eine Familie beschreiben.
Pfam-B liefert automatisch erstellte Cluster, die auf Basis von ClustalW und Hidden
Markov Modellen generiert wurden und familienubergreifend sind [13]. Ausgewiesen
werden dabei nur jene Sequenzen, die mit einem bestimmten Hidden Markov Modell
verglichen eine bestimmte Ahnlichkeit zum Modell aufweisen [13].
10ASTRAL Compendium for Sequence and Structure Analysis - Databases and Tools: http://
astral.berkeley.edu/ (Stand: 02. Juli 2010)11PFAM (Protein Families) Database (Version 24.0, October 2009): http://pfam.sanger.ac.uk/
(Stand: 02. Juli 2010)
3
Grundlagen des Sequenzalignments
Ein Alignment von Sequenzen zielt darauf ab, die in den Sequenzen vorkommenden
DNA-Muster zu identifizieren und deren Position innerhalb der Sequenzen zu bestim-
men. Bei einer Alignierung wird versucht, zwei oder mehrere Zeilen mit je einer Sequenz
- und wenn notwendig durch Einfugen von Leerstellen (Gaps) - so auszurichten, dass es
spaltenweise zu einer bestmoglichen Ubereinstimmung der Symbole kommt. Die Rei-
henfolge der Symbole muss dabei erhalten bleiben.
Abbildung 3.1: Auszug aus einem Alignment von zwei Sequenzen.
Wie aus der Abbildung 3.1 hervorgeht, fuhrt ein Alignment der Sequenzen nicht zwangs-
weise zu einer exakten Ausrichtung identer Symbole. Das Beispiel zeigt nur einige ex-
akte Ubereinstimmungen (orange markiert). Die Berechnung eines Sequenzalignments
entspricht vielmehr einer Suche nach biologisch sinnvollen Kombinationen, die Ahnlich-
keiten in den 3D-Strukturen der Molekule wiedergeben. Erschwerend kommt hinzu,
dass nicht immer eindeutig definierbar ist, nach welchen Regeln einzelne Symbole kom-
biniert werden konnen und wie sinnvoll ein Alignment letztendlich ist. Oft ist das Er-
gebnis nur manuell bewertbar und dient deshalb nur einer Vorverarbeitung fur weitere
Untersuchungen.
17
3. Grundlagen des Sequenzalignments 18
Aufgrund der schnell wachsenden Menge in einschlagigen Datenbanken beschriebener
Sequenzen kommen - trotz der erwahnten Schwierigkeiten - fast nur noch automatisierte
oder computergestutzte Alignmentmethoden zum Einsatz. Beim paarweisen Alignment
kommen sogenannte Substitutionsmatrizen zur Anwendung, welche die Ahnlichkeiten
bekannter Sequenzen abbilden. Beim multiplen Alignment bedient man sich beschrei-
benden Modellen, deren Grundlage maschinell oder (halb-)manuell voralignierte Se-
quenzen und/oder eine vorgeschaltete statistische Auswertung einer großeren Sequenz-
menge ist. Fur eine automatische Ausrichtung zweier oder mehrerer Zeichenketten gibt
es unterschiedliche Verfahren, die sich sowohl in der Empfindlichkeit, als auch in der
Komplexitat der Algorithmen und damit in der Performance deutlich unterscheiden.
Welche Methoden aber auch immer zur Anwendung kommen, alle haben das eine Ziel:
die genetischen Mutationen zu erkennen und die Sequenzen so zu alignieren, dass die
Verwandtschaftsverhaltnisse trotz aller Variationen der Sequenzen sicher erkannt und
die evolutionare Entwicklungen moglichst exakt nachvollzogen werden konnen.
Bevor man sich dem automatischen Alignment widmet, stellen sich einige grundlegende
Fragen, die es prinzipiell und auch im Rahmen dieser Arbeit zu beantworten gilt [21, 32]:
• Wie konnen Sequenzen verglichen werden, auch wenn ihre Symbolabfolgen vorerst
unterschiedlich erscheinen?
• Wie konnen Sequenzabweichungen oder -distanzen bewertet werden?
• Wie kann man Sequenzen so untereinander schreiben, dass die evolutionaren
Vorgange moglichst wiedergegeben werden?
• Wie mussen die Sequenzteile verschoben werden, so dass die untereinander ge-
schriebenen Teile moglichst gut ubereinstimmen?
• Welche Methoden kommen dabei zur Anwendung?
Dieses Kapitel widmet sich der Beantwortung dieser Fragen und fuhrt dabei in die
allgemeinen Grundlagen des Sequenzvergleichs und -alignments ein. Die zwei nachfol-
genden Teile widmen sich zwei speziellen Methoden: dem multiplen Sequenzalignment
und ClustalW und dem Alignment mit Hidden Markov Modellen.
3. Grundlagen des Sequenzalignments 19
3.1 Sequenzvergleich
Die Anordnung und Positionierung der Einzelbausteine, also die Primarstruktur be-
stimmter Proteine, konnen im Laufe der Evolution großen evolutionaren Anderungen
unterworfen sein. Bestimmte Teile werden dabei mehr, andere wiederum weniger kon-
serviert. Die raumlichen und funktionalen Strukturen sind vergleichsweise stabil in
ihrer Entwicklung [21], wohingegen die Sequenzen stark variieren konnen. So konnen
die Aminosauresequenzen innerhalb einer Proteinfamilie aufgrund evolutionarer Ent-
wicklungen betrachtlich voneinander abweichen, obwohl die wesentlichen Eigenschaften
gut erhalten geblieben sind. Die 3D-Strukturen zweier Proteine weisen namlich selbst
dann noch Gemeinsamkeiten auf, wenn die Identitat der Sequenzen mitunter nur noch
30% betragt [32].
Der Vergleich von Sequenzen nimmt eine besondere Stellung ein, weil aus der Ami-
nosauresequenz eines Proteins auf die 3D-Struktur und daruber hinaus auf die Ei-
genschaften eines Proteins geschlossen werden soll. Basis der aus dem Sequenzver-
gleich abgeleiteten Methoden ist die Erkenntnis, dass bei Proteinen eine hohe Se-
quenzahnlichkeit auch eine ahnliche Funktion und/oder 3D-Struktur des Proteins im-
pliziert. Wie kann nun eine Sequenz mit einer anderen vergleichen werden?
Denkt man an einen Sequenzvergleich, so denkt man schnell an einen einfachen, zei-
chenweisen Vergleich zweier Strings. Dabei werden jeweils an der Position n der Strings
A und B die Zeichen an und bn verglichen. Die Sequenzen gelten dann als gleich, wenn
alle Zeichen uber die gesamte Lange der beiden Strings ubereinstimmen. Wie konnen
nun aber Sequenzen miteinander verglichen werden, die einer genetischen Mutation
unterworfen waren, also bei einem zeichenweisen Vergleich deutliche Unahnlichkeiten
aufweisen wurden?
Beim Sequenzvergleich interessiert oft weniger die exakte Folge der Einzelbausteine
der betrachteten Gene, als vielmehr die Interpretation bestimmter Abfolgen und die
Funktionen, die durch eine bestimmte Zusammensetzung der Sequenzen bestimmt wer-
den. Weiß man beispielsweise, dass die Funktionen mehrerer Sequenzen gleicher Fa-
milien ahnlich sind, obwohl die Sequenzen aufgrund ihrer Mutation unterschiedlich
3. Grundlagen des Sequenzalignments 20
sind, so versucht man Muster und Gesetzmaßigkeiten zu erkennen und die Unter-
schiede von Sequenzen zueinander zu quantifizieren. In der Bioinformatik kommen
dafur vielfach statistische Techniken zum Einsatz [15]. Eine Moglichkeit ist beispiels-
weise die Bestimmung des Vorkommens bzw. der Haufigkeit fS(ai) der einzelnen Ami-
nosauren ai oder Codone aus den zu untersuchenden Sequenzen S. Erstellt man daraus
Haufigkeitstabellen, so kann der Unterschied zweier Sequenzen A und B durch die Sum-
me der Haufigkeitsdifferenzen der 20 Aminosauren oder von 61 Sinn-Codone1 errechnet
werden (siehe dazu auch Kapitel 3.2):
Diffa(A,B) :=20∑i=1
|fA(ai)− fB(ai)| (3.1)
Bei der ausschließlichen Betrachtung der Haufigkeiten bleibt die Abfolge der Bausteine
aber unberucksichtigt. Aus diesem Grund ist man auf Methoden ausgewichen, die so-
wohl die Komposition der DNA-Strange als auch die Haufigkeit der Gene berucksichti-
gen und daruber hinaus die Berechnung eines Ahnlichkeitsmaßes moglich machen. Die
in dieser Arbeit verwendeten Hidden Markov Modelle (siehe Kapitel 5) bieten eine
Moglichkeit, Sequenzen zu vergleichen und Ahnlichkeiten numerisch in Form von Ahn-
lichkeitsmaßen auszudrucken [32, 15].
3.2 Ahnlichkeit und Distanz
Vergleicht man zwei Zeichenketten, so kann man mit einem Blick feststellen, ob die-
se gleich oder ungleich sind. Bei gleichen Zeichenketten stimmen alle Zeichen beider
Strings exakt uberein. Weicht nur eines der Zeichen ab oder fehlt ein einzelnes Zei-
chen, so wird man sofort deren Ungleichheit erkennen konnen. Vergleicht man nun
zwei Aminosauresequenzen, die in der Bioinformatik letztendlich auch nur als einfache
Zeichenketten dargestellt werden, so wird ein Vergleich der zu einem Ja/Nein-Ergebnis
fuhrt nicht brauchbar sein. Hier stellt man sich vielmehr die Frage, wie ahnlich zwei
1Bei dieser Methode werden von den 64 moglichen Codone nur 61 sogenannte Sinn-Codone aus-gewertet, die anderen 3 Codone bleiben unberucksichtigt. Diese werden deshalb nicht gezahlt, da siekeinen
”Sinn“ machen; man nennt sie deshalb auch Nonsens-Codon oder auch
”Unsinn“-Codon [25]
3. Grundlagen des Sequenzalignments 21
Sequenzen sind. Die Ahnlichkeit der Zeichen muss also in irgendeiner Form als nume-
rischer Wert ausgedruckt werden konnen.
3.2.1 Hamming-Abstand und -Ahnlichkeit
Nun lage es nahe, bei zwei gleich langen Zeichenketten einfach die Anzahl der uberein-
stimmenden Zeichen zu zahlen. Der aus der Informationstheorie stammende Hamming-
Abstand bedient sich dieser Idee [32]. Der Hamming-Abstand der beiden Sequenzen
PUPPENKISTEN und SUPPENKASPER ergibt 4 und entspricht der Anzahl der vier nicht
ubereinstimmenden Zeichen:
A : P U P P E N K I S T E N| | | | | | | |
B : S U P P E N K A S P E R
Die acht ubereinstimmenden und untereinanderstehenden Zeichen der Sequenzen A
und B sind mit einem ’|’ gekennzeichnet. Um die Hamming-Ahnlichkeit AH zweier
Zeichenketten A := a1a2..an und B := b1b2..bn auszudrucken, wird der Hamming-
Abstand DH(A,B) in Relation zur Gesamtlange n der beiden Zeichenketten gesetzt
und dieser Wert von 1 subtrahiert [23]:
AH(A,B) := 1− DH(A,B)
n(3.2)
Wurden alle Zeichen ubereinstimmen, so entsprache dies dem Wert 1. Kame es zu
keiner Ubereinstimmung so entsprache dies einer Hamming-Ahnlichkeit von 0. Fur die
Sequenzen A und B aus dem obigen Beispiel ergabe die Berechnung der Hamming-
Ahnlichkeit 1− 4/12, also das Ergebnis 0.66.
Die gewichtete Hamming-Ahnlichkeit als erweiterte Variante, berucksichtigt zusatzlich
die einzelnen Positionen innerhalb der Sequenzen. Damit konnen bestimmte Merkmale
einer Sequenz, wie beispielsweise welche Zeichen an welchen Positionen ubereinstimmen
oder nicht ubereinstimmen, uber Merkmalsgewichte in die Berechnung einfließen. Eine
3. Grundlagen des Sequenzalignments 22
Funktion WH(ai, bi, i) ermittelt dabei jeweils das Gewicht der jeweiligen Sequenzzeichen
in Abhangigkeit ihrer Position und Symbolpaarung innerhalb der Sequenzen [23].
AHW(A,B) := 1− DH(A,B)∑ni=1WH(ai, bi, i)
(3.3)
Trotz dieser Verbesserung in Bezug auf Sequenzen erweist sich die Hamming-Ahnlichkeit
beim Sequenzvergleich als nur beschrankt einsetzbar, denn sie kann nur dann berechnet
werden, wenn die beiden Zeichenketten gleich lang sind. Eine Einschrankung, die sie
gerade bei Proteinsequenzen und deren Alignment nahezu unbrauchbar macht. Eine
Kennzahl, bei deren Berechnung die Zeichenketten unterschiedlich lang sein konnen,
ist die sogenannte Levenshtein-Distanz.
3.2.2 Levenshtein-Distanz
Die Levenshtein-Distanz oder der Editierabstand, wie die Levenshtein-Distanz auch ge-
nannt wird, gibt die minimale Anzahl der Editierschritte oder der Kosten an, mit der
eine Zeichenkette in eine andere Zeichenkette uberfuhrt werden kann [30]. Die zu ver-
gleichenden Zeichenketten mussen nicht die gleiche Lange aufweisen. Dies bedingt aber
auch, dass Editieroperationen erlaubt sind, welche eine der beiden Zeichenketten bei
Bedarf verlangert, um die Sequenzen auf die passende Lange und auf Ubereinstimmung
zu bringen.
Das folgende aus [17] entlehnte und leicht abgewandelte Beispiel zeigt die beiden Zei-
chenketten VINTNER und WRITERS und gibt mogliche Operationen an, welche die Zei-
chenkette A zur Zeichenkette B wandeln.
A : V - I N T N E R -| | | |
Op : r i m d m d m m i| | | |
B : W R I - T - E R S
Die Kleinbuchstaben in der Zeile Op symbolisieren dabei die Editieroperationen, welche
das Zeichen an der jeweiligen Position der Zeichenkette A in ein Zeichen der Zeichenket-
te B uberfuhren. Das Zeichen r steht dabei fur die Ersetze- (replace), i fur die Einfuge-
3. Grundlagen des Sequenzalignments 23
(insert), m die Ubereinstimmungs- (match) und d fur die Losch-Operation (delete).
Weist man nun jeder dieser Operationen bestimmte Kosten zu, dann lassen sich daraus
ahnlich wie bei der gewichteten Hamming-Distanz die sogenannte Levenshtein-Distanz
ermitteln [32].
Die Berechnung des Editierabstands der Zeichenketten A := a1a2 . . . am mit der Lange
m und B := b1b2 . . . bn mit der Lange n erfolgt in einer Matrix D mit der Dimension
(m+ 1)× (n+ 1). Beginnend von D0,0 werden die Zellen aufsteigend belegt [30]:
m = |A| n = |B|
∀ 1 ≤ i ≤ m, 1 ≤ j ≤ n
D0,0 = 0; Di,0 = i; D0,j = j
Di,j = min
Di−1,j−1+1ai 6=bj
Di−1,j+1
Di,j−1+1
(3.4)
Nach der Berechnung steht die minimale Anzahl der Editieroperationen in der Matrix-
zelle Dm,n. Die Tabelle 3.1 im Folgeabschnitt auf Seite 26 zeigt eine vollstandig besetzte
Levenshtein-Matrix fur die beiden Zeichenketten VINTNER und WRITERS.
Bei genauer Betrachtung der Levenshtein-Beispiele stellt man fest, dass eine Berech-
nung auch immer eine Alignierung der Sequenzen zur Folge hat. In den beiden Beispiel-
sequenzen A und B wurden hierfur Langenkorrekturen der Zeichenketten vorgenom-
men. Das Zeichen ’-’ symbolisiert dabei eine eingefugte Leerstelle. Symbolpaarungen
konnen so verbessert werden, da die nachfolgenden Symbole eines Strings damit jeweils
um eine Position gegenuber den Symbolen des zweiten Strings verschoben werden. Bei
Versuchen mit einfachen Zeichenketten wird auch schnell deutlich, dass damit durch-
aus unterschiedliche Alignierungen und Levenshtein-Distanzen moglich sind. Weitere
Alignierungen fur A und B waren beispielsweise:
A : V I N T N E R - A : - V I N T N E R -| | | | | | |
Op : r r r m d m m i Op : i r m d m d m m i| | | | | | |
B : W R I T - E R S B : W R I - T - E R S
3. Grundlagen des Sequenzalignments 24
Die beiden Alternativlosungen sind unterschiedlich verlangert und damit auch aligniert
worden. Eine Losung zeigt weniger Lucken auf Kosten mehrerer Fehlpaarungen (Mis-
matches), eine andere wiederum mehr Ubereinstimmungen (Matches), jedoch auf Ko-
sten mehrerer Lucken. Wie ebenso deutlich wird, beeinflussen die Einfugungen (In-
sertions) an unterschiedlichen Stellen die weiteren Editierschritte der Restsequenzen.
Bewertet man nun beispielsweise einen Match mit 0, einen Mismatch mit 1 und eine
Leerstelle mit 2, so kann die Art der Alignierung damit deutlich beeinflusst werden.
Einige Alignierverfahren parametrisieren so ihre Algorithmen.
Aus obiger Erklarung und den Beispielen wird klar, wie eine gewichtete Distanz zweier
Sequenzen berechnet werden kann. Ungeklart bleibt zuweilen noch, wo und wie viele
Lucken im optimalen Fall eingefugt werden sollen. Sucht man ein Verfahren zur Berech-
nung eines optimalen Alignments, so sucht man einen Algorithmus zur Bestimmung
des minimalen Editieraufwands bzw. der minimalen Kosten. Dass ein Probierverfah-
ren hierbei nicht vielversprechend ist, zeigt die Untersuchung der Komplexitat dieser
Aufgabe.
3.2.3 Komplexitat vs. Optimum
Ausgangspunkt der Komplexitatsuntersuchung ist die Festlegung, dass die Ausrichtung
der Symbole und die Positionierung der Lucken dann optimal sind, wenn die Distanz des
dadurch erzeugten Alignments ein Minimum ist. Die einfachste Methode, welche mit
Sicherheit auch zur Optimallosung fuhren wurde ware jene, bei der alle Moglichkeiten
der Alignierung verglichen wurden, um daraus Dmin zu ermitteln.
Schon die genauere Betrachtung eines paarweisen Alignments zeigt die stark sequenz-
langenabhangige Komplexitat eines Brute-Force-Verfahrens (Probierverfahren). Die in
den vorangegangenen Abschnitten vorgestellten Algorithmen versuchen durch geschick-
tes Einfugen von Lucken moglichst viele Paarungen ubereinstimmender oder ahnlicher
Symbole zu erreichen. Im schlechtesten Fall, wenn zwei gleich lange Sequenzen keine
Ubereinstimmungen ermoglichen, muss in den Sequenzen an jeder zweiten Stelle ein
Gap eingefugt werden. Die alignierten Sequenzen wurden sich dadurch in ihrer Lange
addieren [21]. Warum spielt die Sequenzlange dabei eine derart große Rolle?
3. Grundlagen des Sequenzalignments 25
Ein Probierverfahren, bei dem alle moglichen Alignierungen gepruft werden um daraus
die beste Losung zu ermitteln erscheint aussichtslos, wenn man zur Kenntnis nimmt,
wie die Anzahl mGA der moglichen Gap-Anordnungen bei einem paarweisen Vergleich
mit Sequenzen der Lange n kalkuliert wird [21]:
mGA =
(2n
n
)=
(2n)!
(n!)2(3.5)
Wie Hutt und Dehnert in [21] auf Basis dieser Formel berechnen, ergeben sich beim
Vergleich zweier Sequenzen und einer Sequenzlange von n=10 schon 250 Moglichkeiten,
bei n=100 aber schon 1029 mogliche Gap-Anordnungen und zu prufende Falle. Die
Berechnung der minimalen Distanz aus allen moglichen Gap-Anordnungen fuhrt also
schnell an die Grenzen des Moglichen.
Es scheint, als herrsche ein Widerspruch zwischen geringem Aufwand und optimaler
Losung. Dabei handelt es sich bei der am Beginn dieses Abschnitts festgelegten For-
derung um ein klassisches Optimierungsproblem. Vielfach kommt dafur die Technik
des dynamischen Programmierens zum Einsatz. Diese Technik macht zumindest ein
vollstandiges Probierverfahren nicht notwendig, wie das folgende Kapitel zeigt.
3.3 Dynamisches Programmieren
Die dynamische Programmierung kommt vorzugsweise dann zum Einsatz, wenn die
Losung eines algorithmischen Problems uber die Losung von Teilproblemen beschrie-
ben werden kann. Dieses Prinzip von Teile und Herrsche dient als Grundlage fur eine
Reihe von Algorithmen, die sich der Losung eines umfangreichen Problems widmen,
indem sie es in kleinere Teilprobleme zerlegen, die unabhangig voneinander gelost wer-
den konnen. Typische Anwendungen waren jene Spezialfalle, die rekursiv beschrieben
werden konnen. Oft werden aber die Teilergebnisse in Tabellen, Baumen und Listen
zwischengespeichert und spater fur die Losung großerer Probleme wiederverwendet.
Dies erfordert mehr Speicherplatz, ist im Allgemeinen aber effizienter als rekursive
Losungen [41].
3. Grundlagen des Sequenzalignments 26
Dynamisches Programmieren ist auch bei Optimierungsproblemen und im Sequenzali-
gnment weit verbreitet. Es kann zur Losung von Optimierungsproblemen dann erfolg-
reich eingesetzt werden, wenn eine optimale Losung des Problems sich aus optimalen
Losungen der Teilprobleme zusammensetzt. Methoden zur Losung von moglichst guten
Routenplanungen in Navigationssystemen oder in Routern gelten als typische Beispie-
le [41]. Sucht man nun fur zwei Zeichenketten nicht nur irgendeine Alignierung zu
einer Levenshtein-Distanz, sondern jene, bei denen die Operationen in Bezug auf die
Kosten ein Optimum darstellen, so kann ebenso auf die dynamische Programmierung
zuruckgegriffen werden, wie das folgende Beispiel zeigt.
Wendet man die Rechenvorschrift 3.4 fur die beiden Sequenzen VINTNER (A) und
WRITERS (B) an, so kann damit die Matrix D wie in Tabelle 3.1 dargestellt mit den
Symboldistanzen besetzt werden.
Di,j W R I T E R S
0 1 1 2 2 3 3 4 4 5 5 6 6 7 7
V1
1
1 2
2 1
2 3
2 2
3 4
3 3
4 5
4 4
5 6
5 5
6 7
6 6
7 8
7 7
I2
2
2 2
3 2
2 3
3 2
2 4
3 2
4 5
3 3
5 6
4 4
6 7
5 5
7 8
6 6
N3
3
3 3
4 3
3 3
4 3
3 3
4 3
3 4
4 3
4 5
4 4
5 6
5 5
6 7
6 6
T4
4
4 4
5 4
4 4
5 4
4 4
5 4
3 4
5 3
4 5
4 4
5 6
5 5
6 7
6 6
N5
5
5 5
6 5
5 5
6 5
5 5
6 5
5 4
6 4
4 5
5 4
5 6
6 5
6 7
6 6
E6
6
6 6
7 6
6 6
7 6
6 6
7 6
6 5
7 5
4 5
6 4
5 6
5 5
6 7
6 6
R7
7
7 7
8 7
6 7
8 6
7 7
7 7
7 6
8 6
6 5
7 5
4 6
6 4
6 7
5 5
Tabelle 3.1: Matrix zur Berechnung des Levenshtein-Abstands fur die beiden Zeichen-ketten WRITERS und VINTNER. In den Zellen werden die Zwischenergebnisse in klei-ner Schrift dargestellt. Das Minimum der drei Zwischenergebnisse an einer bestimmtenPosition wird jeweils in Normalschrift dargestellt. Der Levenshtein-Abstand kann ausder Zelle Dm,n der letzten beiden Zeichen der Sequenzen abgelesen werden.
3. Grundlagen des Sequenzalignments 27
Wie aus den Levenshtein-Formeln hervorgeht, wird jedes Zeichen einer Sequenz mit
jedem Zeichen der anderen Sequenz verglichen. In den Matrixzellen der Tabelle 3.1
werden jeweils 2 × 2 Ergebnisse der jeweiligen Symbolvergleiche dargestellt. Das Zwi-
schenergebnis oben links in der Zelle errechnet sich aus der Formel fur Di−1,j−1, jenes
rechts oben aus Di−1,j und jenes im linken unteren Bereich der Zelle aus Di,j−1. Zuletzt
wird das Vergleichsergebnis einer Einzelposition rechts unten aus dem Minimum der
drei Zwischenergebnisse ermittelt. Der Levenshtein-Abstand der beiden Zeichenketten
ist das Minimum der summierten Symboldistanzen der letzten beiden Zeichen am und
bn. VINTNER und WRITERS haben demnach, wie aus der Zelle Dm,n abzulesen ist, den
Levenshtein-Abstand 5.
Die Komplexitat zum Befullen der Matrix verglichen mit der Komplexitat eines primi-
tiven Brute-Force-Verfahrens (siehe Gleichung 3.5) zeigt eine deutliche Verbesserung.
Lasst man die Randelemente der Matrix außer acht, so mussen insgesamt m × n Zel-
len berechnet werden, wobei pro Zelle drei Zwischenergebnisse (siehe Gleichungen 3.4)
gerechnet werden, von denen jeweils das Minimum den Zellenwert bestimmt. Demnach
kann die Komplexitat mit O(3mn) beschrieben werden.
Nun dient die gezeigte Matrix nicht nur zur Berechnung der Distanz, sondern auch der
Ermittlung des optimalen Editieraufwands, also jenen mit den minimalen Kosten und
damit ublicherweise der Bestimmung eines optimalen Alignments [32]. Dazu bedient
man sich des sogenannten Ruckwartspfads (traceback path). Verfolgt man beginnend
ab Zelle Dm,n jeweils jenen Weg, durch den das Minimum der einzelnen Zellen be-
stimmt wurde, so lasst sich durch Zuruckverfolgen (backtracking) der Zellen der Pfad
zur Startzelle D0,0 ermitteln. Bestimmen in einer Zelle mehrere Zwischenergebnisse das
Minimum, so entspricht dies einer Verzweigung an dieser Position.
Wendet man dieses Verfahren in der Matrix der Tabelle 3.1 an, so ergeben sich daraus
zwei mogliche Pfade bzw. Alignments:
A : V - I N T N E R - A : V I N T N E R -| | | | | | |
Op : r i m d m d m m i Op : r r r m d m m i| | | | | | |
B : W R I - T - E R S B : W R I T - E R S
3. Grundlagen des Sequenzalignments 28
Legt man nun die Kosten der Operationen mit Matches = 0, Mismatches = 1 und
Inserts = 1 zugrunde und addiert fur die beiden Losungen wieder die Operationsko-
sten, so ergeben sich fur beide Alignments wie erwartet Gesamtkosten in Hohe von
5. Das im Beispiel gezeigte Alignment links benotigt vier Gaps bei vier Matches, das
Alignment rechts nur zwei Gaps bei jedoch nur drei Matches. Die Distanzen beider
Paarungen sind gleich hoch. Beide Alignments sind auf Basis der Kostenfunktionen
gleichwertig und stellen eine optimale Losung fur dieses Alignmentproblem dar.
Die bisher gezeigten Alignments auf Basis der gewichteten Distanzen basieren auf ei-
ner stark vereinfachten Kostenannahme. Matches, Mismatches und das Einfugen von
Lucken werden bei der Bewertung der Distanzen und Alignments mit konstanten Ko-
sten belegt. Die Kosten sind ausschließlich von den lokalen Operationen an deren Sym-
bolpositionen abhangig und vollig unabhangig davon, welche Symbole wie bzw. in wel-
chem Teilabschnitt der Sequenz aligniert werden. Dabei zeigen die Beobachtungen der
Evolutionsprozesse, dass die Vorgange in der Natur weitaus komplexer sind, als in den
obigen Alignmentbeispielen abgebildet. Mit welchen Methoden und Bewertungsverfah-
ren man sich den evolutionaren Vorgaben nahern kann, wird anhand des paarweisen
Sequenzalignments im folgenden Abschnitt eingefuhrt.
3.4 Paarweises Sequenzalignment
Ein paarweises Alignment dient dem Vergleich zweier Sequenzen, um funktionelle und
strukturelle Ahnlichkeiten oder evolutionare Eigenschaften zu erkennen. Aufgrund der
notwendigen Nachbildung evolutionarer Entwicklungen konnen in den zu vergleichen-
den Sequenzen Einfugungen (Insertions), Loschungen (Deletions) oder ein Wechsel
(Mutation) von Symbolen notwendig werden. Fur alle Algorithmen des biologischen
Sequenzvergleichs gilt, dass sie in der Lage sein mussen, diese evolutionaren Vorgange
vollstandig abzubilden und die Alignierungen bewerten zu konnen.
Die in den folgenden Abschnitten vorgestellten Verfahren zur Berechnung globaler und
lokaler Alignments uber die Distanz zweier Zeichenketten zeigen je eine ausgewahlte
Methode, Sequenzen paarweise zu alignieren. Dabei kommt nicht nur der Bewertung
3. Grundlagen des Sequenzalignments 29
der Lucken, sondern auch der Bewertung der Ubereinstimmungen eine besondere Auf-
merksamkeit zu.
3.4.1 Globales Alignment
Ein globales Alignment zweier Sequenzen berucksichtigt alle Symbole der Sequenzen.
Sie kommen vorzugsweise dann zur Anwendung, wenn die zu alignierenden Sequenzen
eine ahnliche Lange aufweisen und große Sequenzubereinstimmungen zu erwarten sind.
Mit einem globalen Alignment versucht man beispielsweise evolutionare Entwicklungen
von Proteinfamilien nachzuvollziehen. Diese eignen sich aufgrund ihrer Homologien gut
fur globale Alignmentmethoden.
Die Berechnung des optimalen globalen Alignments entspricht einem Optimierungs-
problem, welches wie im Abschnitt 3.3 gezeigt uber die dynamische Programmierung
gelost wird. Der dafur verwendete Needleman-Wunsch-Algorithmus [35] ist jenem der
Levenshtein-Distanz (siehe Gleichung 3.4) sehr ahnlich. Anstatt der Konstanten +1
werden den Zellen entsprechend der jeweiligen Operationen symbol- und positions-
abhangige Werte aufaddiert. Die Funktion w(ai, bi) gibt dabei einen von den beiden
Symbolen ai und bi abhangigen Wert (Similarity Score) zuruck, der umso hoher ist, je
ahnlicher die beiden Symbole sind. Die Funktion g ermittelt sogenannte Gap-Strafpunk-
te (Gap-Penalty), also negative Werte, fur das Einfugen einer Lucke an der jeweiligen
Position. Die Art der Berechnung weicht von jener der Levenshtein-Distanz nur insofern
ab, dass pro Zelle anstelle des Minimums jeweils das Maximum der Zwischenergebnisse
weiterverwendet wird. Die in [35] ausschließlich textuelle Beschreibung des Algorithmus
von Needleman-Wunsch lasst sich in wenige Gleichungen zusammenfassen. Die Matrix
M kann wieder von links oben nach rechts unten wie folgt belegt werden [17, 31]:
m = |A| n = |B|
∀ 1 ≤ i ≤ m, 1 ≤ j ≤ n
M0,0 = 0 Mi,0 = i · g(′−′) M0,j = j · g(′−′)
Mi,j = max
Mi−1,j−1+w(ai, bj) : Match (Similarity Score)
Mi−1,j+g(ai,′−′) : Deletion
Mi,j−1+g(′−′, bj) : Insertion
(3.6)
3. Grundlagen des Sequenzalignments 30
Problematisch ist ein globales Alignment dann, wenn zwei Sequenzen A und B mit
deutlich unterschiedlicher Lange aligniert werden sollen. Geht man davon aus, dass
das jeweils gesuchte Symbol der kurzen Sequenz A aufgrund der Beschranktheit des
Alphabets sehr wahrscheinlich an einer der nachstgelegenen Positionen in der langen
Sequenz B gefunden wird und die Positionen in A bis zum Match mit Lucken aufgefullt
werden, dann wird die kurze Sequenz A sich uber den Bereich der langen Sequenz B
”verschmieren“, ohne ein sinnvolles Alignment zu ergeben. Eine negative Bewertung
der Lucken kann dem entgegenwirken, macht aber eine sinnvolle Alignierung trotzdem
kaum moglich [32].
3.4.2 Free-Shift Alignment
Eine erweiterte Form des globalen Alignments ist das Free-Shift Alignment. Dabei wird
eine Folge von Deletions und Insertions am Beginn und am Ende eines Alignments bei
der Bewertung des Alignments nicht berucksichtigt. Diese Lucken am Beginn oder Ende
eines Alignments konnen vermehrt dann auftreten, wenn sich die beiden Sequenzen in
deren Lange deutlich unterscheiden. Verwendet wird diese Alignmentform dann, wenn
uberstehende Prafixes und Suffixes im Ergebnis keine praktische Relevanz haben und
damit unberucksichtigt bleiben konnen.
3.4.3 Lokales Alignment
Ein lokales Sequenzalignment versucht funktions- und strukturrelevante Sequenzteile
zu identifizieren und jene Teile mit einer hohen lokalen Ahnlichkeit zu alignieren. Die
Erweiterung des Needleman-Wunsch- zum Smith-Waterman-Algorithmus ist denkbar
einfach, um diesen fur lokale Alignments anwenden zu konnen [32]. Die erste Zeile
und erste Spalte ist mit 0 zu initialisieren. Zusatzlich ist fur eine Maximierung der
Bewertungen zu sorgen, indem die Gleichung so abgeandert wird, dass die Bewertung
einer Zelle niemals unter Null fallen kann.
3. Grundlagen des Sequenzalignments 31
Nach Anpassung der Gleichungen an obige Forderungen lauten diese [43]:
m = |A| n = |B|
∀ 1 ≤ i ≤ m, 1 ≤ j ≤ n
M0,0 = 0 Mi,0 = 0 M0,j = 0
Mi,j = max
0 : Limitierung auf großer gleich 0
Mi−1,j−1+w(ai, bj) : Match (Similarity Score)
Mi−1,j+g(ai,′−′) : Deletion
Mi,j−1+g(′−′, bj) : Insertion
(3.7)
Hauptsachlich kommen lokale Proteinalignments bei strukturell und funktionell ahn-
lichen Einheiten zum Einsatz, um gemeinsame Domanen von Proteinsequenzen zu iden-
tifizieren [17, 43].
Obwohl sich die Algorithmen zum Generieren eines lokalen und globalen Alignments
unterscheiden, kann keine grundsatzliche und scharfe Grenze in der Anwendung der Me-
thoden und den zu erwarteten Alignments gezogen werden. So ist es durchaus moglich,
dass der fur lokale Alignments verwendete Smith-Waterman Algorithmus aufgrund der
Bewertung der Matches, Mismatches und Lucken ein globales Alignment hervorbringt,
obwohl er besonders bei lokalen Alignments seine Anwendung findet. Ebenso kann ein
Programm, welches den Needleman-Wunsch Algorithmus implementiert, aufgrund der
Bewertung der Lucken oder Parametrisierung auch lokale Alignments generieren. Letzt-
endlich entscheidet nicht nur der Algorithmus, sondern vor allem auch die Bewertung
der Lucken, der Matches und der Mismatches uber das Ergebnis der Alignierung [33].
Schon die Beispiele aus dem Abschnitt 3.3 haben gezeigt, dass durch das Hinzufugen
von Lucken unterschiedliche Alignments moglich sind. Gerade bei Proteinsequenzen
kommt der Bewertung der Leerstellen eine besondere Bedeutung zu [21] und erfordert
deshalb eine genauere Betrachtung.
3. Grundlagen des Sequenzalignments 32
3.5 Lucken und deren Bewertung
Punktmutationen, die zu falschen Paarungen fuhren, durfen je nach Mutation aligniert
werden oder werden durch die Variation der Stringlange durch Lucken (Gaps) aus-
geglichen. Gaps gelten dabei im Allgemeinen als Folge von Einfugungen (Insertions)
oder Loschungen (Deletions) im Laufe der Evolution. Grundsatzlich gilt, je unahnlicher
die relevanten Residuen sind oder je mehr diese uber die Sequenzen verteilt sind, de-
sto mehr Gaps mussen eingefugt werden, um eine Ubereinstimmung der Sequenzen zu
erreichen [38].
Wie aus den Gleichungen 3.6 und 3.7 hervorgeht, fließt die Bewertung einer Leerstelle
jeweils punktuell in die Berechnung ein. Dabei kann eine Leerstelle eine hohe oder eine
niedrige Bewertung haben und die Distanz, vor allem aber das Alignment damit deut-
lich beeinflussen. Wird eine Lucke schlechter bewertet als ungleiche Symbolpaarungen,
so wird der Paarung der Vorzug gegeben. Kommt den Lucken eine bessere Bewertung
zu als etwaigen Falschpaarungen, so werden an diesen Positionen vermehrt Lucken
eingefugt werden, um Mismatches zu vermeiden. Das Einfuhren oder Fortsetzen einer
Lucke wird mit Strafpunkten (Gap-Penalties) bzw. negativen Werten verrechnet. Eine
einfache Bewertungsregel konnte beispielsweise lauten:
Match = 3 Mismatch = 0 Gap-Penalty = -1
Nun zeigt sich aber gerade bei Proteinen, dass diese in mehreren kurzeren Teilsequen-
zen hohe Ahnlichkeiten aufweisen, da sie haufig aus identischen Domanen zusammen-
gesetzt sind. Damit kann es erforderlich werden, an jenen Stellen in der Sequenz, an de-
nen durch evolutionare Mutationen zusatzliche Domanen eingefugt wurden, bei einem
Alignment entsprechend langere Lucken einzufugen [32]. Eine konstante Bewertung
der Gap-Penalties kann diese Vorlage der Natur nicht abbilden. Um dieser Forderung
gerecht zu werden, ist eine diesen Entwicklungen angepasste Bewertung der Lucken ein-
zufuhren. Moglich ist dies durch affine Kostenfunktionen, welche die Bindungsstarke
der Proteinsequenzen abbilden.
Bei der sogenannten linearen Bewertung eines Gaps der Lange n, fließt uber die Glei-
chung f(n) = −n·PGap die Starke der Gap-Penalty PGap in die Kostenfunktion f(n) ein
3. Grundlagen des Sequenzalignments 33
[21]. Bei dieser Bewertung wird aber nicht zwischen dem Einfuhren und dem Fortset-
zen einer Lucke unterschieden. Anders ist dies bei der Bewertung der Lucken mithilfe
affiner Gap-Penalties. Hierbei wird das Einfuhren einer Lucke mit dem Gap-opening
Penalty und das Verlangern einer Lucke um eine Leerstelle mit dem Gap-extension
Penalty bewertet [42]. Die Gesamtkosten einer Lucke errechnen sich dann aus den
beiden gesondert betrachteten Gap-Penalties und der Lange der Lucke wie folgt [32]:
f(n) := PGapOpening + n ·GGapExtension (3.8)
In der Regel unterscheiden sich diese Penalties deutlich und das Gap-opening Penalty
wird schlechter, also mit einem hoheren jedoch negativen Penalty bewertet, als das
Gap-extension Penalty. Ublicherweise werden beim Alignment von Proteinsequenzen
Werte von -5 bis -19 fur das Einfuhren einer Lucke veranschlagt und -1 bis -3 fur das
Verlangern einer Lucke um eine Position [32]. Ist die Differenz der beiden Penalties groß,
das Offnen einer Lucke also”teuer“ und die Verlangerung einer Lucke vergleichsweise
”gunstig“, so werden wenige, stattdessen aber langere Lucken eingefugt, als wurden die
Penalties etwa gleich bewertet werden.
3.6 Substitutionsmatrizen
Die in den vorangegangen Abschnitten vorgestellten Algorithmen zur Berechnung der
Distanz zweier Sequenzen und eines globalen oder lokalen Alignments arbeiten auf Ba-
sis von Scores, welche sowohl einen Match als auch einen Mismatch und das Einfugen
einer Lucke bewerten. Die Bewertung von Matches oder Mismatches wird auf Basis
statistischer Auswertungen getroffen, bei der aus einer geeigneten Menge von Protei-
nen die relativen Haufigkeiten von Aminosauren und deren Substitutionshaufigkeiten
bestimmt werden. Diese Haufigkeiten werden in Scores transformiert und damit Sub-
stituationsmatrizen2 gefullt [21].
2Substitutionsmatrizen werden manchmal auch (similarity matrix ) oder Austauschmatrix (muta-tion data matrix ) genannt, das sie die Ahnlichkeit von Aminosauren innerhalb einer Sequenz und dasMutationsverhalten bewerten.
3. Grundlagen des Sequenzalignments 34
Eine Substitutionsmatrix beschreibt mit einer Menge von Scores sa,b die Wahrschein-
lichkeiten der Ersetzung der Aminosaure a durch eine andere Aminosaure b (und um-
gekehrt) in einer bestimmten Sequenz. Dabei wird die Ubereinstimmung mit einer
selten vorkommenden Aminosaure hoher bewertet als die Ubereinstimmung mit einer
haufig vorkommenden Aminosaure. Zudem versucht man aufgrund der Eigenschaften
der Aminosauren, den Mismatch zweier funktionell ahnlicher Aminosauren hoher zu
bewerten als den Mismatch zweier funktionell verschiedener Aminosauren. [32]
Substitutionsmatrizen speichern fur jedes Paar bzw. alle moglichen Aminosaurekombi-
nationen nicht deren Wahrscheinlichkeiten, sondern deren davon abgeleitete Scores. Um
dies zu verdeutlichen wird manchmal auch die Bezeichnung Scoring-Matrizen vorgezo-
gen. Die Zellen der Hauptdiagonale beschreiben jeweils das Zusammentreffen zwei-
er gleicher Aminosauren (Match), alle anderen Felder eine Fehlpaarung (Mismatch)
und die Substitution von einer Aminosaure zu einer anderen. Die einfachsten Scoring-
Matrizen sind Identitatsmatrizen, in denen alle Diagonalelemente, jene die einen Match
anzeigen, den gleichen Score enthalten und alle anderen Elemente den gleichen und
niedrigeren Score enthalten. Damit wird zum Ausdruck gebracht, dass die Wahrschein-
lichkeit einer Mutation zu einer beliebigen anderen Aminosaure konstant ist [32, 21].
Eine Vereinfachung, die in der Praxis kaum noch zur Anwendung kommt.
Eine Betrachtung der Gleichungen 3.6 und 3.7 im Zusammenhang mit der Bewertung
von Matches und Mismatches uber Substitutionsmatrizen macht deutlich, dass diese
zusammen mit der Bewertung von Gaps einen wesentlichen Einfluss auf die Art des
Alignments haben. Liegen die Scores der Matches und Mismatches dicht beisammen, so
wird ein kompakteres Alignment gebildet. Ist ein Match mit einem deutlich positiveren
Score bewertet als ein Mismatch, das heißt ein Match kompensiert mehrere Mismatches,
so werden die Alignments in Lange gezogen und schwacher ausgepragt sein, sofern das
Einfugen von Lucken nicht ausdrucklich negativ bewertet wird (siehe Abschnitt 3.5).
Ziel aller Scoring-Matrizen ist eine moglichst reale Abbildung der Evolution. Dies ist
auch der Grund, dass fur die jeweilige Anwendung eine passende Matrix ausgewahlt
werden muss. Zwei beim Alignment von Aminosauren haufig verwendete Matrizen sind
die sog. PAM-Matrizen und die BLOSUM-Matrizen.
3. Grundlagen des Sequenzalignments 35
3.6.1 PAM-Matrizen
Die 1978 von Dayhoff et al. in [10] eingefuhrten PAM-Matrizen (PAM steht fur Point
Accepted Mutations) basieren auf sogenannten PAM-Einheiten, mit denen die evolu-
tionare Distanz zweier Aminosauresequenzen gemessen wird. Zwei Sequenzen S1 und
S2 unterscheiden sich per Definition um eine PAM-Einheit, wenn S2 wahrend seiner
evolutionaren Entwicklung durch eine Reihe akzeptierter Punktmutationen aus S1 ent-
standen ist und dabei pro 100 Residuen durchschnittlich eine Punktmutation auftrat.
Unter einer akzeptierten Punktmutation versteht man in diesem Zusammenhang den
Austausch einer Aminosaure durch eine andere und die Weitervererbung dieser Mutati-
on. Berucksichtigt werden nur jene Mutationen, welche die Funktion des Proteins nicht
oder zu seinem Vorteil verandern. Jene Mutationen, die durch Insertion oder Deletion
entstanden sind, bleiben unberucksichtigt. Aufgrund der Definition wurde man vermu-
ten, dass Sequenzen, die durch 100 PAM-Einheiten divergieren, zu 100 % voneinander
abweichen. Aufgrund wechselnder Mutationsereignisse werden im Laufe der Evolution
jedoch immer wieder Veranderungen vorgenommen, so dass zuvor eingefuhrte Mutatio-
nen zu einem spateren Zeitpunkt durchaus wieder ruckgangig gemacht werden konnen.
Aus diesem Grund sind PAM-Werte auch uber 100 moglich. Sequenzen mit 250 PAM-
Einheiten konnen durchaus noch Ubereinstimmungen von 20 % und mehr aufweisen
[11, 32].
PAM-Matrizen quantifizieren fur alle Aminosauren die Wahrscheinlichkeiten einer evo-
lutionaren Entwicklung, die zu bleibenden Mutationen fuhren. Die in Abbildung 3.2
auszugsweise und im Anhang A.1 vollstandig dargestellte PAM250-Matrix zeigt die
Wahrscheinlichkeiten einer evolutionaren Distanz von 250 PAM-Einheiten.
Tabelle 3.2: Auszug aus einer PAM250 Matrix
Um die Werte verstandlicher zu machen, wurden die Elemente mit 100 multipliziert,
damit handelt es sich bei den Werten in der Tabelle um Prozentangaben. Die Tabelle
3. Grundlagen des Sequenzalignments 36
zeigt beispielsweise, dass bei einer Aminosauresequenz S1, an deren erster Position ein
A (Ala) stand, mit einer Wahrscheinlichkeit von 13 % auch in der mutierten Sequenz
S2 ein A (Ala) stehen wird. Ebenso besteht eine 3%ige Chance, dass in der mutierten
Sequenz S2 ein A (Ala) zu einem R (Arg) mutiert.
Bis heute werden PAM-Matrizen zum Alignieren von Sequenzen verwendet, obwohl die
Frage nach der evolutionaren Distanz von Proteinen bis heute nicht klar beantwortet
werden kann. Welche PAM-Matrix am ehesten anzuwenden ist, hangt auch von der
jeweiligen Fragestellung und von den Ausgangssequenzen ab. Eine pragmatische Vor-
gehensweise ist dabei die Verwendung mehrerer unterschiedlicher Matrizen. Matrizen,
die sich vor allem bei der Identifikation von Proteindomanen bewahrt haben, gehoren
zur Familie der BLOSUM-Matrizen [32].
3.6.2 BLOSUM-Matrizen
Die spater als die PAM-Matrizen entwickelten BLOSUM-Matrizen (BLOSUM steht da-
bei fur BLOck SUbstitution Matrix ), gehen auf die Untersuchungen von Henikoff und
Henikoff [19] zuruck. Die Autoren verwendeten fur ihre Untersuchungen Daten aus
einer BLOCKS Datenbank. Dabei werden stark konservierte Regionen aus multiplen
Sequenzalignments extrahiert und in die Datenbank aufgenommen, wobei Lucken dabei
nicht zugelassen sind. Die Funktion der Motive muss nicht bekannt sein. Die Auswer-
tung der BLOSUM-Blocke erfolgt spaltenweise, indem fur jede Aminosaure a ∈ A und
b ∈ B erst deren Haufigkeit f(a) und f(b) und anschließend die Haufigkeit f(a, b) der
Aminosaurenpaarungen a und b in samtlichen Spalten berechnet wird. Die Berechnung
der Scores erfolgt als Log-odds-Verhaltnis in der Form [32]:
sa,b := logsf(a, b)
f(a)f(b)(3.9)
Die Eintrage in der Odd-Score-Matrix reprasentieren damit Wahrscheinlichkeiten. So-
mit kann der Score eines Gesamtalignments uber die Einzelwahrscheinlichkeiten be-
rechnet werden. Die Verwendung des Logarithmus bietet den Vorteil, dass die Werte
miteinander addiert werden konnen und nicht multipliziert werden mussen, wie bei
3. Grundlagen des Sequenzalignments 37
Wahrscheinlichkeiten ublich. Auf das Endergebnis hat diese Transformation keinen
Einfluss. Die Basis s des Logarithmus ist in der Regel so gewahlt, dass der Score nur
kleine Werte einnimmt, sich meist im Bereich von -10 bis +10 bewegt und eine Run-
dung der Werte auf Integer keinen großen Einfluss auf das Endergebnis hat. Oft wird
die Basis 2 benutzt.
Fuhrt man dieses Verfahren wie in Gleichung 3.9 beschrieben durch, so kann man
nachvollziehen, dass die Auswahl der Sequenzblocke einen wesentlichen Einfluss auf die
Große des Scores hat. Werden uberproportional viele einander sehr ahnliche Blocke ver-
wendet, so werden diese das Ergebnis weitgehend dominieren und seltene Mutation nur
einen geringen Einfluss auf die Scores haben. Eine Verbesserung des Verfahrens besteht
nun darin, einander sehr ahnliche Sequenzblocke beim Scoring bewusst auszuschließen.
Dadurch bekommen die Mutationen der unahnlicheren Sequenzen ein großeres Gewicht,
was sich in einem hoheren Score fur seltenere Paarungen niederschlagt. Diese Verbes-
serung des Modells wird in der Bezeichung der BLOSUM-Matrizen dokumentiert. Eine
BLOSUM62-Matrix wie in Tabelle 3.3 dargestellt, wird nur aus Sequenzblocken berech-
net, deren Sequenzen beim paarweisen Vergleich maximal 62 % identische Residuen
aufweisen. Ublich sind daruber hinaus noch BLOSUM50- oder BLOSUM80-Matrizen.
Die Verwendung der BLOSUM-Matrix zeigt das folgende einfache Beispiel. Gegeben
sind die beiden Sequenzblocke A mit der Sequenz LVLHVWAK und B mit der Sequenz
LCMKSLEH, die den zwei alignierten Sequenzen d1a6ma_a.1.1.2 und d1asha_a.1.1.2
entnommen wurden:
A: Leu Val Leu His Val Trp Ala LysL V L H V W A K| | | | | | | |L C M K S L E H
B: Leu Cys Met Lys Ser Leu Glu His
Der Score S des Alignments der Sequenzen A und B errechnet sich aus der Summe
der Einzelwerte fur jede Symbolpaarung, die der BLOSUM-62-Tabelle 3.3 entnommen
werden konnen:
S: 4 - 1 + 2 - 1 - 2 - 2 - 1 - 1 = -2
3. Grundlagen des Sequenzalignments 38
Tabelle 3.3: BLOSUM-62 Scoringmatrix fur 20 Aminosauren fur das Alignment vonProteinsequenzen.
Wie in den Abschnitten 3.2 zum Thema Ahnlichkeit und Distanz gezeigt wird, kann
uber die Summe der Bewertung einzelner Symbolpaarungen oder der Lucken eine Be-
wertung der Ahnlichkeiten ganzer Sequenzen erfolgen. Diese Ahnlichkeitswerte symbo-
lisieren nicht nur die Distanz zweier Sequenzen, sondern konnen auch zur Konstruktion
eines Alignments herangezogen werden. Bisher wurden diese Werte wie dimensionslose
Großen betrachtet, ohne sie auf ein Bezugssystem zu beziehen. Wie diese Großen in
Bezug auf ein Normal verwendet werden konnen und wie sie zu interpretieren sind,
erklart der folgende Abschnitt.
3.7 Scoring von Alignments
Scoring steht im Zusammenhang zum Sequenzalignment fur ein Bewertungsverfah-
ren, mit dem Ahnlichkeiten quantifiziert werden. Die Ahnlichkeit oder der als Distanz
bezeichnete Abstand einer Sequenz zu einem Bezugssystem wird durch numerische
3. Grundlagen des Sequenzalignments 39
Werte bzw. Scores ausgedruckt. Je hoher der Score einer Sequenz ist, desto großer
ist die Ahnlichkeit der Sequenz zum Bezugssystem. Die Suche nach einem optimalen
Alignment ist de facto immer der Suche nach einem hohen oder maximalen Score gleich-
zusetzen, der eine gute oder im besten Fall optimale Ubereinstimmung von Sequenzen
beschreibt.
Die in den vorangegangenen Kapiteln eingefuhrten Methoden alignieren Sequenzen und
platzieren Lucken so, dass eine moglichst hohe Ubereinstimmung bei Symbolpaarungen
erreicht wird. Die Distanz der Paarungen ist dabei moglichst gering, da die Menge der
ubereinstimmenden oder ahnlichen Paarungen optimiert wird. Bei einem paarweisen
Sequenzvergleich wird zwar letztlich die Distanz zweier Sequenzen ausgedruckt und die
Basiseinheit der Scores wird uber die Scoring-Matrizen bestimmt, das Bezugssystem
wird dabei jedoch nicht festgelegt3. Erst eine Schatzung eines Bezugssystems uber
empirische Methoden kann einen Anhaltspunkt dafur liefern, was einen guten Score
ausmacht.
3.7.1 Empirische Festlegung eines Bezugssystems fur Scores
Ein Score fur sich alleine betrachtet ist wertlos, solange keine Vergleichsmoglichkeiten
und kein Bezugssysytem zur Verfugung stehen. Mit den bisher diskutierten Abstands-
werten ist ein Vergleich von Sequenzabstanden nur untereinander moglich, eine Frage
ist bisher aber unbeantwortet geblieben: Wie hoch muss ein Score sein, um eine signi-
fikante Ahnlichkeit zweier Sequenzen auszudrucken? Um diese Frage beantworten zu
konnen, wird auf zufallige (randomisierte) Sequenzen zuruckgreifen [21].
Ausgangspunkt der Betrachtungen ist die Einsicht, dass bei einem Vergleich zweier
beliebiger Sequenzen auch zufallige Ubereinstimmungen auftreten konnen. Verfahren,
die das Einfugen von Gaps erlauben, vergroßern die Menge der ubereinstimmenden oder
ahnlichen Symbolpaarungen zusatzlich, da die Gaps mit dem Ziel platziert werden, eine
moglichst hohe Ubereinstimmung zu erreichen. Halt man die Symbole einer Sequenz
nun bei und erlaubt gleichzeitig eine zufallige Positionierung der Symbole (Shuffling),
3Uber die wahrscheinlichkeitstheoretische Betrachtung der Log-odds im Bayes Framework lasstsich so ein Bezugssystem prinzipiell beschreiben. Derartige Berechnungen sind im Allgemeinen jedochnumerisch schwer handhabbar.
3. Grundlagen des Sequenzalignments 40
so erhalt man eine randomisierte Sequenz. Wiederholt man diesen Vorgang, so konnen
die generierten Zufallssequenzen einer weiteren Analyse zugefuhrt werden.
Mit der sogenannten Signifikanzanalyse wird das Alignment zweier Sequenzen mit dem
Alignment randomisierter Sequenzen verglichen. In der Datenanalyse und Modellierung
bezeichnet man damit einen Vergleich mit kunstlich generierten Ersatzdaten (Surro-
gatdaten), die als Stellvertreter fur die Originaldaten dienen und diese reprasentativ
ersetzen. Die Frage die damit beantwortet werden soll ist, ob Sequenzen signifikante,
also wesentliche Ahnlichkeiten aufweisen, oder ob die Unterschiede in den Sequenzen
zufalliger Natur sind. Zu diesem Zweck werden aus den Originaldaten kunstliche Da-
tensatze, meist randomisierte Datensatze erzeugt, diese Daten einer Analyse zugefuhrt
und die Ergebnisse dieser Daten mit den Originaldaten verglichen. Entscheidend ist
dabei, inwieweit die Originaldaten von dem durch die Surrogatdaten beschriebenen
Ergebnisbereich entfernt liegen. Angewendet auf den Vergleich randomisierter Sequen-
zen bedeutet dies, dass die erste Sequenz beibehalten wird, die zweite Sequenz jeweils
randomisiert wird, um daraus ein optimales Alignment zu erhalten. Jene randomisierte
Sequenz mit dem maximalen Score entspricht dem optimalen Alignment und beant-
wortet gleichzeitig die Frage, wie groß der optimale (maximale) Score ist [21].
Eine vollstandige Sequenzrandomisierung wie oben beschrieben, ist in der Praxis zu
aufwendig. Die Kenntnis uber die Score-Verteilung lasst hierbei aber einige Vereinfa-
chungen zu.
3.7.2 Verteilung der Scores
Ist die Verteilungsfunktion der Scores bekannt, so konnen mit weit weniger Aufwand
als die Randomisierung der Sequenzen die Parameter der Verteilungsfunktion geschatzt
werden und damit die Signifikanz des ermittelten Scores belegt werden. Karlin und
Altschul haben 1990 in [24] die theoretische Verteilung der Scores zufalliger Sequenzen
hergeleitet und empirisch gezeigt, dass die Verteilungsfunktion der Scores sowohl im
Falle der Alignments mit als auch ohne Gaps annahernd einer Extremwertverteilung
folgt [21].
3. Grundlagen des Sequenzalignments 41
Wiederholt man die oben beschriebene Randomisierung und Alignierung entsprechend
oft mit anderen randomisierten Sequenzen, so erhalt man eine Haufigkeitsverteilung der
Scores auf der Basis der beiden Ausgangssequenzen. Das in Abbildung 3.2 dargestellte
Histogramm beschreibt, welche Scores in welchem Umfang aufgrund der Sequenzzusam-
mensetzungen zu erwarten sind. Ohne hierbei auf die zwei konkreten Ausgangssequen-
zen und deren Scores naher einzugehen, lasst sich aus dem dargestellten Histogramm
eine Verteilungsfunktion erkennen.
Abbildung 3.2: Histogramm der Score-Verteilung randomisierter Sequenzen. Der ro-te Balken zeigt den Wert des ausgewahlten Originalalignments. Die Kurve zeigt dietheoretische Verteilung nach Karlin-Altschul (basierend auf [21])
So zeigt das Histogramm in der Abbildung 3.2, dass der Score von -11 die großte
Haufigkeit hat. Der Score des Originalalignments am außeren rechten Rand der Score-
Verteilungskurve (siehe roter Balken) zeigt eine Ahnlichkeit der beiden Sequenzen an,
die uber die Sequenzzusammensetzung deutlich hinausgeht. Ein Score von 11 weist auf
einen Score hin, dessen zufalliges Auftreten bei der dargestellten Verteilung als unwahr-
scheinlich gelten kann. Die in der Abbildung dargestellte theoretische Verteilungskurve,
lasst das Maximum in Bereich von etwa -8 erwarten.
Theoretisch waren auch andere Verteilungsfunktionen moglich, wenn die Berechnung
der Scores auf anderen Bewertungsgrundlagen als jene der Symbolhaufigkeiten und
deren Paarungen basiert. Die Kenntnis der Verteilungsfunktion der Scores ist aber
notwendig, soll die Große eines Scores einem Bezugssystem zugeordnet werden konnen.
3. Grundlagen des Sequenzalignments 42
3.7.3 Standardisierung der Scores
Aus dem Histogramm geht augenscheinlich hervor, dass mit einem Score von 11 eine
uber dem Zufall liegende Ahnlichkeit der beiden Originalsequenzen wahrscheinlich ist.
Diese Erkenntnis soll jedoch auch quantitativ ausgedruckt werden konnen. Als einfach
zu errechnender Wert gilt dafur der Z-Score, den beispielsweise auch der FASTA4-
Algorithmus verwendet. Der Z-Score wird dabei einer Standardisierung (Z-Transforma-
tion) uber den Mittelwert und der Standardabweichung unterzogen, um die Scores
miteinander vergleichbar zu machen und zu standardisieren [37].
Der Z-Score gibt die Abweichung eines konkreten Wertes s vom Mittelwert S der
Scoreverteilung als Vielfache der Standardabweichung σS an [32].
Z-Score :=s− SσS
(3.10)
Je großer der Score s ist, desto mehr weicht er an den rechten Rand der Verteilung,
also in die positive X-Richtung aus. Umso weiter der Score rechts vom Mittelwert S
liegt, desto großer ist die Wahrscheinlichkeit, dass der Score keinem Zufall entspricht
bzw. das Alignment nicht die Folge einer zufalliger Ubereinstimmungen ist.
Als zum Z-Score alternative Kenngroßen kommen auch der P-Wert (P-Value) und
E-Wert (E-Value) zur Anwendung. Der P-Wert gibt die Wahrscheinlichkeit an, dass
ein gleich hoher oder hoherer Score durch eine zufallige Ubereinstimmung zustande
gekommen ist. Der E-Wert hangt von der Menge der Eingangssequenzen ab. Er gibt
die erwartete Zahl jener Sequenzen der gesamten Menge an, die beim Vergleich mit
einer randomisierten Sequenz zum errechneten oder einen hoheren Score fuhren wurde
[21].
Die bisherigen Betrachtungen haben sich ausschließlich paarweisen Sequenzvergleichen
und -alignments gewidmet. Die multiple Alignierung von Sequenzen als Erweiterung
obiger Verfahren gilt algorithmisch als wesentlich komplexer. Eine genauere Betrach-
tung dieser Problemstellung ist Thema des folgenden Kapitels.
4Der FASTA-Algorithmus wurde 1985 von David J. Lipman und William R. Pearson als FASTPfur Proteine entwickelt [28]; die aktuellen FASTA-Tools sowohl fur Proteine als auch Nukleotide sindunter http://www.ebi.ac.uk/Tools/fasta/index.html zu finden.
4
Multiples Sequenzalignment
Unter einem Multiplen Sequenzalignment (MSA) versteht man die Erweiterung eines
paarweisen Alignments auf drei oder mehrere Sequenzen. Dabei werden die Sequenzen
- aquivalent einem paarweisen Alignment - so zueinander ausgerichtet, dass die Distanz
der einzelnen Sequenzen moglichst gering ist und am Ende alle Sequenzen die gleiche
Lange haben. Notigenfalls mussen dafur Gaps in einzelne Sequenzen eingefugt werden,
die ein lokales Alignment der Sequenzen uberhaupt erst moglich machen [38].
Ein multiples Sequenzalignment A besteht wie ein paarweises Alignment auch aus dem
Alphabet ΣA und um das fur eine Lucke symbolisierende Zeichen”-“. Es gilt wie schon
beim paarweisen Alignment, dass keine Spalte nur aus Lucken bzw. den Zeichen”-“
bestehen darf.
Aus der Forderung des gleichzeitigen Alignierens mehrere Sequenzen und dem dadurch
deutlich steigenden Aufwand gegenuber dem paarweisen Sequenzalignment, entste-
hen zusatzliche algorithmische Herausforderungen, die in diesem Kapitel eine kurze
Einfuhrung finden. Nachdem erlautert wird, wo ein multiples Sequenzalignment seine
Anwendung findet und den Aufwand rechtfertigt, muss - wie beim paarweisen Alignie-
ren schon - eine Methode zur effizienten Suche nach dem optimalen Alignment gefunden
werden.
43
4. Multiples Sequenzalignment 44
4.1 Anwendung multipler Sequenzalignments
Proteine konnen in sogenannte Proteinfamilien zusammengefasst werden. Ausschlag-
gebend dafur ist die Ahnlichkeit der Proteine bezuglich ihrer Proteinstruktur oder
ihrer Funktionalitat. Proteine die aus der gleichen Proteinfamilie stammen, enthalten
dieselben charakteristischen Domanen, was wiederum auf enge evolutionare Zusam-
menhange der Proteine hinweist [32]. Beobachtet man uber viele Sequenzen hinweg
lokale Ahnlichkeiten zwischen den jeweiligen Sequenzen, so handelt es sich dabei oft
um strukturelle Ahnlichkeiten der jeweiligen aus einer Familie stammenden Protei-
ne [21]. Damit man nun strukturell ahnliche Sequenzen und globale Ahnlichkeiten
erkennen kann, muss die Suche und der Vergleich der Sequenzen auf eine großere Se-
quenzmenge ausgedehnt werden [32]. Ein paarweiser Vergleich zweier Sequenzen ist
dafur ungeeignet, da das Gesamtsystem und die Gesamtcharakteristik einer Familie
dabei unberucksichtigt bleiben. Ein multiples Sequenzalignment deckt in vielen Fallen
die wesentlichen Informationen der evolutionaren Entwicklung und der funktionalen
Zusammenhange auf, bei denen ein paarweiser Vergleich versagen wurde [21].
Bei Sequenzen mit geringen Ahnlichkeiten ist die Erkennung der uber die gesamte
Sequenz verteilten relevanten Residuen mit statistischen Methoden schwierig, da die
Signale durch die Sequenzvarianten in einem Rauschen zu verschwinden drohen. Oft-
mals fuhrt nur die Alignierung ganzer Sequenzfamilien zu besseren Ergebnissen, denn
erst der Vergleich mehrerer Sequenzen verstarkt die Signale der konservierten Residuen
so, dass sie vom einfachen Signalrauschen unterschieden werden konnen [32].
Multiple Sequenzalignments bieten wesentliche Vorteile bei der Erkennung und Bewer-
tung konservierter und variabler Positionen innerhalb von Proteinfamilien, doch gerade
das Alignieren mehrerer Sequenzen zu einem MSA gilt als algorithmisch herausfor-
dernd und aufwendig [32]. Einerseits sucht man nach einer optimalen Alignierung aus
einer Unzahl moglicher Sequenzanordnungen, andererseits steigt der Aufwand mit der
zunehmenden Anzahl an Sequenzen uberproportional an. Eine bioinformatische Kern-
frage der letzten Jahre war und ist, wie mehrere Sequenzen mit einem uberschaubaren
Aufwand zu einem multiplen Sequenzalignment (MSA) aligniert werden konnen [21]
4. Multiples Sequenzalignment 45
und wie dabei mit dem Widerspruch optimale Losung zu algorithmischen Aufwand
umgegangen wird.
4.2 Komplexitat multipler Sequenzalignments
Dass ein Probierverfahren nur beschrankt zweckmaßig ist, um eine optimale Alignie-
rung zweier Sequenzen zu erreichen, wurde schon im Kapitel 3.2.3 deutlich. Die Be-
rechnung des maximalen Score aus allen moglichen Gap-Anordnungen fuhrt schnell
an die Grenzen des Moglichen, da die algorithmische Komplexitat und der daraus
generierte Rechenaufwand eines exakten Vergleichs so hoch sind, dass ein derartiges
Vorgehen trotz immer schnellerer Hardware in der Praxis kaum durchfuhrbar ist. In-
folge der Definition der dynamischen Programmierung (siehe Kapitel 3.3) skaliert der
Rechenaufwand mit dem Produkt der beiden Sequenzlangen. Vergleicht man die Aus-
gangssequenzen mit den Sequenzen einer Datenbank, so skaliert der Aufwand linear
mit der Anzahl der zu vergleichenden Sequenzen [21]. Bei einem multiplen Sequenz-
alignment nimmt der Aufwand noch einmal deutlich zu, indem jede Sequenz mit jeder
anderen Sequenz verglichen werden muss.
Betrachtet man den Aufwand exakter iterativer Methoden bis hin zu den schon deut-
lich schnelleren Methoden auf Basis der dynamischen Programmierung, so ist zu erwar-
ten, dass bei einem multiplen Sequenzalignment mit einem direkten Vergleich ahnliche
dynamische Methoden wie bei einem paarweisen Vergleich kaum zur Losung des Pro-
blems herangezogen werden konnen. Hutt und Dehnert [21] quantifizieren die Grenzen
des Machbaren, indem sie die iterativen Methoden des paarweisen Sequenzvergleichs
auf das multiple Sequenzalignment umlegen und den Aufwand diskutieren. Bei N Se-
quenzen mit annahernd gleicher Lange L benotigt man LN Score-Eintrage in einer
N -dimensionalen Matrix um daraus einen optimalen Score zu berechnen. Der uber dy-
namische Programmierung ermittelte Pfad der optimalen Alignmentkonstruktion fuhrt
durch einen N -dimensionalen Raum beschrieben in einer N -dimensionalen Struktur.
Ohne die schrittweise Herleitung der Autoren wiederzugeben1 wird am Endergebnis
1Die Quantifizierung des Aufwands und die Verallgemeinerung des mathematischen Problems samteiner ausfuhrlichen Herleitung findet sich in [21, S. 179-183]
4. Multiples Sequenzalignment 46
doch klar,”dass eine unmittelbare Verallgemeinerung der Algorithmen des paarweisen
Alignments kein gangbarer Weg zur Konstruktion eines multiplen Alignments ist“ [21,
S. 183].
Die Große heutiger Datenbanken und der Aufwand eines paarweisen Vergleichs zweier
Sequenzen legen nahe, die Anzahl der Sequenzen und damit die Sequenzvergleiche
moglichst fruhzeitig zu minimieren.
Grundsatzlich gibt es drei wesentliche Einflussgroßen, die den Rechenaufwand eines
multiplen Sequenzvergleichs bestimmen:
• die Lange der jeweiligen Sequenzen bzw. Regionen
• die Große des Alphabets
• die Anzahl der zu untersuchenden Sequenzen
Diese drei Faktoren konnen vorweg nur schwer oder nicht reduziert werden. Die Lange
der Sequenzen kann nicht beliebig reduziert werden, da diese großtenteils”naturgege-
ben“ sind. Gene haben in der Regel eine Lange von mehr als 80 Codone [32], die es
bei einem globalen Alignment zu alignieren gilt. Gleiches gilt fur die Große des Al-
phabets, welches ebenso nicht reduzierbar ist, da bestimmte Sequenzen naturgemaß
hohere Variationen im Code aufweisen und somit mit vielen unterschiedlichen Codone
beschrieben werden.
Gerade beim multiplen Sequenzalignment spielt die Menge der Ausgangssequenzen ei-
ne wesentliche Rolle. Schon in der Einfuhrung zu Kapitel 4 wird festgestellt, dass beim
multiplen Sequenzalignment die Anzahl der Eingangssequenzen bewusst hoch gehalten
wird, um zu qualitativ besseren Ergebnissen zu gelangen. Eine Methode, die bei einem
multiplen Sequenzalignment deshalb nahezu unumganglich ist, betrifft die Vorauswahl
von Sequenzen aus einer ursprunglich großeren Sequenzmenge. Wenn es gelingt, aus
einer zunachst großen Menge von Sequenzen mit einem effizienten Verfahren eine rele-
vante Teilmenge zu extrahieren, deren Sequenzen wahrscheinlich gut zu alignieren sind
oder umgekehrt, deren Sequenzen man bei den Betrachtungen mit hoher Wahrschein-
lichkeit ausschließen kann, so kann der Aufwand damit deutlich reduziert werden. Dies
wirft jedoch die Fragen auf, wie und welche Sequenzen aus einer Menge vorselektiert
4. Multiples Sequenzalignment 47
werden konnen und wie und in welcher Reihenfolge diese dann zu einem MSA ali-
gniert werden konnen. Iterative bzw. progressive Methoden versuchen dieses Problem
in effizienter Weise zu losen.
4.3 Iterative Generierung eines MSA
Die Komplexitatsbetrachtungen im Abschnitt 4.2 haben gezeigt, dass ein ’Probierver-
fahren’, das aus allen moglichen Losungen die beste Alignierung sucht, zur Bestim-
mung multipler Alignments kaum zielfuhrend sein kann. Eine effizientere Methode ist
die iterative Erstellung eines MSA auf Basis eines paarweisen Vergleichs. Dabei wer-
den aus einer Menge von Sequenzen die Scores aller Paarungen mittels dynamischer
Programmierung berechnet. Am Beginn wird aus der Menge alle Paarungen jene Paa-
rung mit dem hochsten Score aligniert. Dieses Sequenzpaar bildet ein erstes MSA. In
jedem folgenden Iterationsschritt wird eine weitere Sequenz dem bestehenden MSA
hinzugefugt. Ausgewahlt wird dabei immer jene der verbleibenden Sequenzen, die dem
bisher generierten MSA am nachsten ist. Wie diese Distanzabschatzung genau erfolgt,
ist applikationsabhangig. Das im folgenden Abschnitt als Beispiel angefuhrte ClustalW
leitet vom paarweisen Abstand der Sequenzen binare Baume ab, die in einem weiteren
Alignierungsschritt die Reihenfolge festlegen. Im Laufe der Alignierung wird der Baum
dann Schritt fur Schritt abgearbeitet [33, 32].
An das MSA hinzugefugt wird jede neue Sequenz, indem sie zu den ursprunglichen
Sequenzen im MSA, jedoch ohne Berucksichtigung der Lucken, aligniert wird. Muss
aufgrund der Alignierung der neuen Sequenz eine Lucke in eine Sequenz der MSA
eingefugt werden, so muss diese Lucke in alle Sequenzen des MSA eingefugt werden.
Die Iterationsschleife wird solange fortgefuhrt, bis entweder alle Sequenzen zum MSA
hinzugefugt wurden, oder bis die Scores der verbleibenden Sequenzen nicht mehr den
Anforderungen genugen.
ClustalW, einer der meistbenutzten Algorithmen zur Konstruktion eines multiplen Se-
quenzalignments [46] basiert auf einem iterativen Verfahren. Clustal generiert paarweise
Alignments und fugt diese in einem weiteren Arbeitsgang zu einem MSA zusammen
[33]. Wie Clustal im Detail ablauft, beschreibt der folgende Abschnitt.
4. Multiples Sequenzalignment 48
4.4 Multiples Sequenzalignment mit ClustalW
Clustal ist ein progressiver Algorithmus und ein gleichnamiges darauf basierendes Soft-
warepaket zur Berechnung eines multiplen Sequenzalignments. Die Applikation steht
sowohl als Shell-Applikation als ClustalW als auch mit einer grafischen Benutzerober-
flache als ClustalX fur Linux, Mac und Windows zur Verfugung2.
Der grundlegende Clustal Algorithmus gilt, obwohl er in der ersten Fassung schon zu
Beginn der 1990er Jahre beschrieben wurde und 1994 Thompson et al. in [45] die
erweiterte Version unter ClustalW veroffentlicht haben, bis heute als eine der Stan-
dardmethoden zur Erstellung multipler Sequenzalignments [33]. ClustalW ist die Wei-
terentwicklung von Clustal, einem Algorithmus zur paarweisen Alignierung von Se-
quenzen, der hier nicht naher beschrieben wird3. Das W im Namen steht fur gewichtet
(Weighted) und weist darauf hin, dass Sequenzen entsprechend deren Ahnlichkeiten zu
anderen Sequenzen mehr oder weniger gewichtet werden konnen, um die Alignierung
der Sequenzen zu MSAs positiv zu beeinflussen [33, 45].
ClustalW nutzt ein progressives Alignment, um ein MSA schrittweise zu berechnen.
Der grundlegende Multiple Alignment Algorithmus lauft dabei in drei Schritten ab, die
in [45] und [32] wie folgt beschrieben sind:
1. Im ersten Schritt werden fur alle Paarungen separat globale Alignments be-
stimmt. Zudem werden mittels dynamischer Programmierung uber einfache Sco-
res jeweils die Distanzen jeder Paarung berechnet und diese in einer Distanzma-
trix abgelegt. Der Begriff einfach bezieht sich dabei auf den Berechnungsaufwand.
So werden optional Gaps beim Scoring nur mit konstanten Werten berucksichtigt,
um den Aufwand der Alignierung zu reduzieren. Eine Vereinfachung, die vor allem
bei großen Datenmengen genutzt wird.
2ClustalW und ClustalX Multiple Sequence Alignment - Multiple alignment of nucleic acid andprotein sequences. Siehe unter: http://www.clustal.org/
3Auf der Webprasenz http://www.clustal.org/ der Clustal-Entwickler finden sich unter denReferenzen eine ausfuhrliche Zusammenfassung der wesentlichen Veroffentlichen, die auch eine stetigeWeiterentwicklung der Software bis heute dokumentieren.
4. Multiples Sequenzalignment 49
2. Auf Basis der generierten Distanzmatrix erfolgt im nachsten Schritt die Berech-
nung eines phylogenetischen Baumes. Mit Hilfe des Neighbour-Joining-Algorith-
mus wird ein binarer Baum generiert, der einerseits die Reihenfolge der iterativen
Alignierung der einzelnen Sequenzen zu einem MSA festlegt, und andererseits
auch die Gewichtung der einzelnen Sequenzen bzw. deren Score bestimmt. Das
Gewicht der Sequenzen wird dabei uber die Distanz des Endknotens zur Wurzel
(root) des Baums ermittelt.
3. Im letzten Schritt wird die Alignierung auf Basis des erstellten Clusterbaumes
vorgenommen. Die Alignierung wird mit jener Paarung begonnen, welche den
hochsten paarweisen Score aufweist. Diese Paarung ist auch jene, die als erster in
den Baum eingefugt wurde. In weiterer Folge werden die verbleibenden Sequenzen
sukzessiv und in der Reihenfolge des Hinzufugens der Knoten in den Baum dem
Alignment hinzugefugt. Weist ein Knoten zwei Kindsequenzen auf, so werden die-
se vorher aligniert und dann erst dem MSA hinzugefugt. Hat ein Knoten nur eine
Kindsequenz oder ein oder mehrere Alignments, so werden diese wenn notwendig
zuvor zusammengefasst und dem MSA hinzugefugt. Ausdrucklich hervorzuheben
ist, dass die Alignierung an den Endknoten immer mit der Originalsequenzen
erfolgt, und nicht mit den etwaig eingefugten Gaps der vorangegangenen Alignie-
rung zur Bestimmung des binaren Baums.
Der ClustalW Algorithmus ist verglichen mit anderen Methoden effizient, da eine paar-
weise Alignierung nur einmal berechnet werden muss. Dennoch hat ClustalW auch ei-
nige Nachteile. Lokale Sequenzahnlichkeiten bleiben beispielsweise oft unterbewertet,
da der Algorithmus auf eine Berechnung der Distanzen auf Basis eines globalen Ali-
gnments basiert. Sequenzmengen, die nur eine oder wenige lokale Gemeinsamkeiten
besitzen, sich ansonsten aber weitgehend unterscheiden, sind fur ClustalW wenig ge-
eignet.
Die Reihenfolge bei der Alignierung spielt ebenso wie die Distanz der alignierten Se-
quenzen vor allem am Beginn der Iteration eine wesentliche Rolle. Zu stark divergente
Sequenzen gerade zum Beginn der Iteration, als auch die Bewertung der Einfuhrung
von Lucken und der Verlangerung der Lucken beeinflussen das Ergebnis empfindlich.
4. Multiples Sequenzalignment 50
Ist die Ahnlichkeit der ersten alignierten Sequenzen gering, so ist die Wahrscheinlich-
keit hoch, dass es zu falschen Paarungen kommt oder Lucken an falschen Positionen
eingefugt werden. Nachteilig ist dabei besonders, das einmal durchgefuhrte Alignments
nachtraglich nicht mehr korrigiert werden. Feng und Doolittle beschreiben diese Eigen-
art mit”Once a gap, always a gap“ [14].
Gerade bei stark divergenten Sequenzen ist die Wahl der Parameter wesentlich und
beeinflusst die Ergebnisse in hohem Maße [32]. Auf die jeweiligen Ausgangssequen-
zen angepasste Substitutionsmatrizen und eine sinnvolle Bewertung der Gaps konnen
dem Abdriften des Algorithmus entgegenwirken. ClustalW unterstutzt eine Reihe un-
terschiedlicher Protein Weight Matrizen (BLOSUM 30, PAM 350 und einige mehr),
um den Alignierungsprozess an die Sequenzen anzupassen. Fehlalignments und falsche
Ergebnisse durch lokale Minima sind trotzdem nicht auszuschließen und bleiben im
ungunstigen Fall unentdeckt. Der Algorithmus liefert, wie die meisten anderen Algo-
rithmen auch, keinen Gesamtscore, der eine objektive Bewertung der Qualitat des MSA
moglich machen wurde.
4.5 Bewertung eines multiplen Sequenzalignments
Fast alle Verfahren bleiben die Bewertung des multiplen Sequenzalignments durch die
Berechnung eines Gesamtscores schuldig. Damit kann die Qualitat eines MSA weder
verglichen, noch zum Ausdruck gebracht werden. Methoden der Alignmentbewertung,
wie sie bei paarweisen Alignierungen verwendet werden, konnen nicht ohne weiteres
auf ein MSA umgelegt werden. Bei einer geringen Sequenzanzahl k mit einer geringen
Sequenzlange n ist die Berechnung auf Basis der Scores paarweise alignierter Sequenzen
theoretisch moglich. Dabei werden die Scores aller moglichen Sequenzpaarungen eines
MSA aufsummiert. Merkl und Waack weisen in [32] aber ausdrucklich darauf hin, dass
in der Praxis die Berechnung uber die Sum of Pairs wenig zielfuhrend ist, da der Auf-
wand mit O(nk) exponentiell zur Anzahl der Sequenzen uber alle vertraglichen Maße
steigt. Die Abschatzung des Aufwands mit ublichen Großen zeigt schnell die Grenzen
dieser Bewertungsmethode auf. Dies ist auch ein Grund dafur, dass im Abschnitt 7
alternative Bewertungsmethoden zur Anwendung kommen.
5
Hidden Markov Modelle
Untersucht man eine Menge genetisch verwandter Sequenzen, um sowohl die Haufigkeit
ihres Auftretens als auch das Mutationsverhalten der Symbole mit probabilistischen
Modellen zu beschreiben, so stellt man fest, dass die Anordnung der Symbole wie zu
erwarten nicht zufallig ist. Die Ordnung gehorcht bestimmten Regeln, welche letztend-
lich das Muster der evolutionaren Entwicklungen beschreiben. Eine der wesentlichen
Aufgaben der Mustererkennung in der Bioinformatik besteht darin, biologisch zusam-
menhangende Fragmente in einer langeren Sequenz zu erkennen oder Ahnlichkeiten
zwischen verschiedenen Sequenzfragmenten festzustellen. Diese Fragmente konnen als
Signale verstanden werden, deren statistische Eigenschaften es in Modellen zu beschrei-
ben gilt.
An der Vielfalt der genetischen Mutationen ist schon zu erkennen, dass die Regeln
der Evolution keinesfalls trivial und statisch sind. Sie lassen sich nicht mit einfachen
Wenn-dann-Regeln nachbilden, sondern bedurfen komplexen Modellbeschreibungen.
Bei der Beschreibung der Modelle stutzt man sich deshalb auf Wahrscheinlichkeiten, die
eine wahrscheinlich richtige Vorhersage zulassen und so in der Auswahl der Sequenzen
bessere Ergebnisse als reine Zufallstreffer moglich machen.
Eines in der Muster- und vor allem Spracherkennung weitverbreiteten, beschreibenden
Modelle auf Basis von Wahrscheinlichkeiten ist das Hidden Markov Modell (HMM).
Das nach A. Markov1 benannte statistische Modell wurde am Beginn der 1990er von
1Andrei Andrejewitsch Markow, 1856 - 1922, russischer Mathematiker; wesentliche Beitrage zurWahrscheinlichkeitstheorie und Analysis; verfasste die Theorie der stochastischen Prozesse (Markow-Prozesse), die in den Hidden Markov Ketten zur Anwendung kommen [1].
51
5. Hidden Markov Modelle 52
Haussler et al. [18] zur Untersuchung von Proteinsequenzen in die Bioinformatik ein-
gebracht. Das Modell wird durch zwei stochastische Prozesse beschrieben: durch Zu-
standsubergange und Emissionen, die gemeinsam die Markovsche Kette beschreiben.
5.1 Markov-Ketten
Aminosauresequenzen werden in der Bioinformatik ublicherweise als Zeichenketten dar-
gestellt. Vergleicht man die Sequenzen untereinander, so stellt man fest, dass sowohl
die Anordnung der Symbole als auch die Lange der Sequenzen variieren und zufallig
erscheinen. Markovsche Ketten beschreiben die Wahrscheinlichkeiten P(Y = y), mit
denen eine Zufallssequenz y uber dem Alphabet ΣA = {A,C,D, . . . , V,W, Y } eine
beliebige, feste Anordnung der Symbole y=y1y2y3 . . . yn ∈ ΣnA annimmt.
Mit der Festlegung, dass die Lange L ∈ N der Zufallssequenz diskret und ebenso zufallig
ist, ergibt sich eine Wahrscheinlichkeit P(L = n), dass die Sequenz Y eine bestimmte
Lange n hat. Das Produkt der Wahrscheinlichkeit P(L = n) und der Wahrscheinlich-
keit fur das Auftreten einer bestimmten Sequenz y mit der Lange n ergibt die totale
Wahrscheinlichkeit [32]:
P(Y = y, L = n) = P(L = n) · P(Y = y|L = n) (5.1)
Die Berechnung der Wahrscheinlichkeit P(L = n) bedarf in diesem Zusammenhang
keiner weiteren Erlauterungen. Die Ermittlung der Wahrscheinlichkeiten, dass die Zu-
fallssequenz Y mit einer zuvor festgelegten Lange n eine bestimmte Symbolabfolge y
ergibt, entpuppt sich bei genauerer Betrachtung jedoch als in der Praxis zu aufwen-
dig. Die Wahrscheinlichkeit P(Y = y|L = n) errechnet sich, indem der Prozess zur
Erzeugung der Zufallssequenz als eine Folge von stochastischen Einzelprozessen ange-
sehen wird, die Zufallselemente Y1,Y2,. . . generieren. Zu jedem Zeitpunkt n nimmt das
Zufallselement Yn einen Wert yn aus ΣA an. Betrachtet man den Prozess zu einem
spateren Zeitpunkt n, so haben zu diesem Zeitpunkt die Zufallselemente Y1,Y2,. . .,Yn−1
bereits die Werte y1,y2,. . .,yn−1 angenommen. Die Folge dieser Uberlegung ist eine
iterative Anwendung der Einzelwahrscheinlichkeiten, aus denen sich die gemeinsame
5. Hidden Markov Modelle 53
Wahrscheinlichkeit der Zufallselemente Y1,Y2,. . .,Yn ermitteln lasst[32]:
P(Y = y|L = n) := P(Y1 = y1, Y2 = y2, . . . Yn = yn) (5.2)
Unter der Verwendung bedingter Wahrscheinlichkeiten wird der iterative Charakter
der Gleichung noch deutlicher.
P(Y1 = y1,Y2 = y2, . . . Yn = yn) =
P(Y1 = y1)·P(Y2 = y2|Y1 = y1) · · ·
· · ·P(Yn = yn|Y1 = y1, Y2 = y2, . . . , Yn−1 = yn−1)
(5.3)
Dieser iterative Ansatz in Form der Zerlegung auf Einzelereignisse gilt fur jedes Wahr-
scheinlichkeitsmodell einer Sequenz [21]. Die Gleichung zeigt, dass mit zunehmender
Anzahl n die Berechnung der Wahrscheinlichkeiten immer aufwendiger wird und die
Anzahl der Moglichkeiten exponentiell wachst. In der Praxis ist diese Berechnungsme-
thode demnach kaum anwendbar.
Markovsche Ketten begrenzen per Definition die Anzahl der Auswahlmoglichkeiten.
Voraussetzung dieser Begrenzung ist die Annahme, dass eine eingeschrankte Kenntnis
der vorangegangenen Symbolanordnung ebenso gute Prognosen uber die nachfolgenden
Symbole zulasst wie bei Kenntnis aller der aktuellen Position vorangegangenen Sym-
bolanordnungen des aktuellen Prozesses. Eine Einschrankung des Modells auf n · |ΣA|
Auswahlmoglichkeiten bedeutet, dass die Auswahl der Symbole fur das Zufallselement
Yn zum Zeitpunkt n nur vom Zeitpunkt selbst und vom Wert yn−1 ∈ ΣA abhangt, den
das vorgehende Zufallselement Yn−1 angenommen hat [21, 32]:
P(Yn = yn|Y1 = y1, Y2 = y2, . . . , Yn−1 = yn−1) = P(Yn = yn|Yn−1 = yn−1) (5.4)
Man bezeichnet dieses Modell der Gleichung 5.4 als Markovsche Kette erster Ordnung2.
Die bedingte Wahrscheinlichkeit erstreckt sich dann nur noch uber ein Symbol und
nicht uber die gesamte vorangegangene Sequenz.
2Markov-Ketten hoherer Ordnung haben eine großere Reichweite. Deren bedingte Wahrscheinlich-keiten erstrecken sich dann uber mehrere Symbole.
5. Hidden Markov Modelle 54
Markov-Ketten konnen als endliche Automaten interpretiert werden, deren Zustande
mit gerichteten Kanten verbunden werden, die mit Ubergangswahrscheinlichkeiten ver-
sehen sind [21].
Abbildung 5.1: Erweiterter Zustandsgraph einer Markov-Kette der Menge Σ = {(a, l) :a ∈ A, l = 1, 2, . . . , n} mit der Lange n = 4 und dem Startknoten (aus [32])
Verfolgt man im obigen Beispiel der Abbildung 5.1 die gerichteten Kanten der Markov-
schen Kette vom Startknoten bis zum Zeitpunkt n, so wird deutlich, dass der erweiter-
te Zustandsgraph der von einem Zustand zum nachsten fuhrt, einer zufalligen Irrfahrt
gleichkommt. Wesentlich ist die Festlegung, dass die Knoten nur in der vorgegebenen
Richtung der Kanten erreicht werden konnen. Die Knotenmenge besteht aus der Menge
{Start} ∪ {ΣA × . . .× ΣA}. Am Beginn der Irrfahrt befindet man sich im Startknoten
und fahrt nun zufallig, jedoch der Wahrscheinlichkeitsverteilungen entsprechend den
gerichteten Kanten entlang von einem Knoten yi zum nachsten yi+1. Die Schatzung
der Wahrscheinlichkeit, die Knoten in einer bestimmten Reihenfolge anzufahren, er-
gibt sich aus dem Produkt des sich durch die Irrfahrt ergebenden Zustandsgraphs.
Die erweiterte Form einer Markov-Kette findet in der Bioinformatik in Hidden Markov
Modellen ihre Anwendung.
5.2 Hidden Markov Modelle
Ein Hidden Markov Modell (HMM) beschreibt einen stochastischen Prozess, der uber
zwei Zufallsprozesse definiert wird [21, 32]:
• Der erste stochastische Prozess beschreibt die internen Zustandsubergange zwi-
schen beliebigen und diskreten Zustanden aus dem internen Zustandsraum Λ =
5. Hidden Markov Modelle 55
{1, 2, . . . , λ} auf Basis der Ubergangswahrscheinlichkeiten von einem Zustand in
den anderen. Er entspricht einer Markov-Kette. Dieser Prozess beschreibt die in-
terne Zustandsfolge x=x1,x2,. . .,xl des Modells. Sie ist nicht beobachtbar bzw. ist
”versteckt“ (engl. hidden), daher auch die Bezeichnung Hidden Markov Modell.
Eine Folge von Zustanden x∈ Λl mit der Lange l und dem Startzustand x0 wird
als Pfad bezeichnet.
• Der zweite Prozess generiert entsprechend einer zustandsabhangigen Wahrschein-
lichkeitsverteilung zu jedem Zeitpunkt i eine sichtbare Emission y=y1,y2,. . .,yl
aus dem Emissionsalphabet Σ = {1, 2, . . . , σ}. Die Folge der Emissionen ist be-
obachtbar und entspricht beispielsweise bezogen auf Sequenzen den einzelnen
Aminosauren. Eine Folge von Emissionen y∈ Σl mit der Lange l wird Beobach-
tung genannt3.
Der Hidden Markov Prozess kann - gleich der Markov-Kette - als zufallige Irrfahrt
in einem erweiterten Zustandsgraphen des Modelles aufgefasst werden. Bei einem ge-
dachten Experiment stehen in jedem Knoten entlang des Zustandsgraphen zwei Au-
wahlmoglichkeiten zur Verfugung. Die ersten Auswahlmoglichkeit bildet die Emissionen
aus dem Emissionsalphabet entsprechend den jeweiligen Emissionswahrscheinlichkeiten
ab. Die zweite stellt die gerichteten Kanten zu den Nachbarknoten zur Auswahl, wobei
hierbei nur die Knoten des Zeitpunkts i+ 1 erreicht werden konnen.
5.3 Profil-HMMs
In der Bioinformatik kommen vorwiegend Profil-HMMs zur Anwendung. Der Aufbau
des Profil-HMMs leitet sich von den Sequenzen ab, mit denen das Hidden Markov Mo-
dell trainiert wird. Beim Training des Modells werden die Sequenzen eines MSA spal-
tenweise analysiert und die Verteilung der Symbole innerhalb der Spalten berechnet.
Spalten mit großen Symbolahnlichkeiten bilden einen Konsens (engl. Consensus) fur
die im multiplen Alignment enthalten Proteinteile. Da in einigen Sequenzen Teile feh-
len konnen oder zusatzlich zu den Consensus-Spalten Teilfolgen in einzelne Sequenzen
eingefugt wurden (siehe Kapitel 3), bestehen multiple Sequenzalignments neben den
3In der konkreten Anwendung mit Aminosauresequenzen ware Σ = ΣA
5. Hidden Markov Modelle 56
Konsenspositionen auch aus Einfugungen (Insertions) und Loschungen (Deletions), die
im HMM berucksichtigt werden mussen.
Profil-HMMs werden mit drei grundsatzlichen Zustandstypen beschrieben [32]:
• Match-Zustande mi (in der Abbildung 5.2 mit Quadraten symbolisiert) beschrei-
ben eine Position i in der Sequenz, die innerhalb des Modells zum Konsens gehort.
Dies muss nicht zwingend bedeuten, dass zwei exakt gleiche Symbole der Menge
Σ ubereinstimmen. Es konnen durchaus mehrere und unterschiedliche Symbole
zum Konsens gehoren, wenn die aus einem MSA einer Proteinfamilie bestehenden
Trainingsdaten (Profil-MSA) die Haufung mehrerer unterschiedlicher Symbole an
der betreffenden Position bestatigen. Befindet sich der Prozess im Zustand mi,
so emittiert er mit einer der Spalte i zugrundeliegenden Wahrscheinlichkeit ein
Symbol aus dem Alphabet ΣA.
• Insert-Zustande ii (in der Abbildung 5.2 als Rhomben dargestellt) erzeugen ein
zusatzliches Symbol in der Sequenz, obwohl mit den Modellspalten kein Kon-
sens erreicht wurde. Der Prozess emittiert ein Symbol aus dem Alphabet ΣA.
Nachdem die Trainingsdaten selbst keine Emissionswahrscheinlichkeiten fur die
betreffende Position beschreiben, kann beispielsweise auf Hintergrundwahrschein-
lichkeiten zuruckgegriffen werden. Voraussetzung dafur ist, dass fur jedes Symbol
der Menge ΣA eine Wahrscheinlichkeit großer Null definiert ist.
• Bei Delete-Zustanden di (in die Abbildung 5.2 als Kreise eingebunden) wird im-
mer eine Lucke, also der Buchstabe”-“ emittiert. In diesem Fall gilt die Spalte
als verarbeitet.
Um das Modell zu vervollstandigen, wird es um den initialen Zustand und den ter-
minalen Zustand erweitert. Der initiale Zustand 0 entspricht dem Startzustand, der
zu einem spateren Zeitpunkt niemals mehr angenommen werden kann. Der termina-
le Zustand ∞ emittiert nichts und ist immer der letzte Zustand der erreicht werden
kann. Der Prozess beginnt immer im initialen Zustand und endet immer im terminalen
Zustand.
5. Hidden Markov Modelle 57
Abbildung 5.2: Schematische Darstellung eines Profil-HMM mit der Lange n (aus [32])
Die gerichteten Kanten in der Abbildung 5.2 symbolisieren die jeweiligen Zustandsuber-
gange (Transitionen) von einem Zustand in einen anderen. Sie zeigen auch an, dass sich
der Zustand entlang einer Kante nur in eine Richtung verandern kann. Jede dieser Tran-
sitionen ist mit positionsabhangigen Ubergangswahrscheinlichkeiten verknupft. Diese
geben an, mit welchen Wahrscheinlichkeiten der Pfad den Weg uber eine bestimmte
Kante nimmt. Zusatzlich wird in den Insert- und Delete-Zustanden gemaß der Emissi-
onsverteilung ein Symbol aus der Menge ΣA emittiert.
In diesem Zusammenhang sollten auch mogliche Probleme bei der Parametrisierung
eines Profil-HMMs mit einem Profil-MSA nicht unerwahnt bleiben [32]:
• Es ist durchaus problematisch, wenn aus einer dunnbesetzten Spalte der Trai-
ningsdaten Wahrscheinlichkeiten geschatzt werden sollen. Diese Spalten sind mit-
unter wenig reprasentativ.
• Die Schatzung der Ubergangswahrscheinlichkeiten zu und aus Insertionszustanden
kann sich dann als problematisch erweisen, wenn die Trainingsdaten kaum Inser-
tionszustande aufweisen.
Fur beide Probleme besteht die Losung darin, Hintergrundwahrscheinlichkeiten in die
Schatzung einfließen zu lassen, darauf basierende Pseudocounts einzufuhren die verhin-
dern, dass eine Wahrscheinlichkeit Null wird, oder die Spalte selbst als Ergebnis einer
Einfugeoperation zu interpretieren und nicht in das HMM aufzunehmen [12].
5. Hidden Markov Modelle 58
5.4 Verwendung von Profil-HMMs
In obigen Abschnitten wurden mehrmals die Parallelen eines Hidden Markov Modells
mit einem Automaten, der nach bestimmten Regeln das Modell in Einzelschritten
durchlauft, deutlich gemacht. In jedem dieser Schritte, die sich aus der Konfiguration
des Modells ergeben, wird ein Symbol emittiert. Welcher Pfad durch das Modell ge-
nommen wird und welche Symbole emittiert werden, bestimmen sowohl die Ubergangs-
als auch die Emissionswahrscheinlichkeiten des Modells.
Bisher ist unerwahnt geblieben, welche Fragestellungen mit einem HMM beantwortet
werden konnen und wie ein Profil-HMM fur ein Sequenzalignment konkret genutzt
werden kann. Der Umgang mit dem Profil-HMM lasst sich in drei Aufgabenbereiche
untergliedern [21]:
• Training : Das Training des Modells erfolgt uber die Schatzung der Modellpa-
rameter mit einer Menge alignierter Sequenzen. Aus der Zusammensetzung der
Sequenzen und deren Alignierungen ergibt sich letztlich die Modellbeschreibung.
• Evaluation: Sie entspricht der Gewinnung einer Sequenz aus dem internen Pfad.
• Decodierung : Mit der Decodierung wird der interne Pfad einer gegebenen (beob-
achteten) Sequenz ermittelt.
Im Rahmen dieser Arbeit gilt den Aufgaben Training, vor allem aber der Decodierung
das besondere Interesse. Anhand eines konkreten Beispiels werden beide Aufgaben im
Folgenden naher beschrieben, um am Ende des Kapitels zu zeigen, wie der Zustandspfad
zur Alignierung einer Sequenz mit einem MSA fuhrt.
5.4.1 Training eines Profil-HMM
Das Training eines Profil-HMM entspricht der Schatzung der HMM-Parameter auf
Basis von Sequenzen mit einer bekannten internen Struktur. Die Sequenzen werden
aufgrund deren Strukturahnlichkeiten zu einem MSA aligniert und bilden so die Da-
tengrundlage des Trainingsprozesses. Wahrend des Trainingsvorgangs werden in ei-
nem iterativen Prozess die Emissions- und Transitionswahrscheinlichkeiten unter dem
5. Hidden Markov Modelle 59
geschatzten Modell maximiert [21] (siehe auch Maximum Likelihood Schatzer in [12]
und [32]).
Wie werden Sequenzen eines Eingangs-MSA interpretiert, um uberhaupt die Eingangs-
daten fur eine Schatzung zu erhalten? Krogh et al. [26] entwickelten dafur ein Losungs-
verfahren, das wie folgt zusammengefasst werden kann [12, 32]:
Gegeben ist ein multiples Sequenzalignment Am×l = (aij) mit m Sequenzen und l
Spalten. Sie enthalten Buchstaben aus ΣA und Gaps”-“ aus vorangegangenen Alignie-
rungen. Um das Profil-HMM mit dem gegebenen MSA A zu trainieren, wird aus den
einzelnen Sequenzen des MSAs jeweils ein Pfad konstruiert.
Dazu wird jenen Spalten ein Match-Zustand M zugeordnet, bei der die Anzahl der
Buchstaben aus ΣA großer einem bestimmten Schwellenwert ρ ist. Die Anzahl der
Match-Zustande uber alle Spalten des MSAs bestimmt die Lange n des Modells MA,
bestehend jeweils aus einem Match- mi, Insertion- ii und Deletion-Zustand di. Spalten
mit entsprechender Besetzung werden in das Modell aufgenommen und zu Modellspal-
ten (in dieser Arbeit mit der Markierung M versehen). Dunnbesetzte Spalten werden
als Insertions interpretiert und nicht in das Modell aufgenommen (diese Spalten werden
hier mit einem N markiert). Das Beispiel in Abbildung 5.3 zeigt die Zuordnung der
Spalten als eine MN-Reihe uber den Modellsequenzen.
Mit dem HMM kann fur jede Sequenz ai,· im MSA der Pfad mit der Lange L ermittelt
werden. Ist eine Spalte einem Match-Zustand M zugeordnet, so gilt das Zeichen aik der
Sequenz in der jeweilige Spalte k als vom Match-Zustand mi emittiert, vorausgesetzt
es handelt sich nicht um das Zeichen”-“ fur ein Gap; andernfalls wird das Zeichen
dem Deletion-Zustand di zugeschrieben. Ist das Zeichen aik keinem Match-Zustand M
zugeordnet, so gilt das Zeichen vom Insertion-Zustand ii emittiert, wenn es zuvor schon
eine Spalte gab, die einem Match-Zustand m zugeordnet wurde; andernfalls gilt das
Zeichen als von i0 emittiert.
5. Hidden Markov Modelle 60
5.4.2 Decodierung eines Profil-HMM
Zu jeder beobachteten Sequenz kann aus dem Profil-HMM durch Decodierung die ge-
meinsame Wahrscheinlichkeit eines Pfades aus ΣHMM und der Emission aus ΣA berech-
net werden. Je hoher diese Wahrscheinlichkeit ist, als desto geringer gilt die Distanz
der Beobachtung zum Modell. Gleiches gilt fur Scores, denn HMM-Scores werden in
den meisten Fallen von den Modell-Wahrscheinlichkeiten abgeleitet [32].
Aus der schematischen Darstellung des Profil-HMMs in Abbildung 5.2 wird deutlich,
dass es aufgrund der nicht beschrankten Lange der durch das HMM modellierten Se-
quenzen unendlich viele und demnach unendlich lange Pfade im Modell gibt, die zu
einer Emission gleich einer beobachteten Sequenz fuhren. Die Bildung eines optima-
len Alignments fuhrt uber die Suche nach dem wahrscheinlichsten aller Pfade. Die
Methode, bei der alle Pfade berechnet werden, um jenen Pfad mit der hochsten ge-
meinsamen Wahrscheinlichkeit zu erhalten, erscheint impraktikabel, da die Zahl der
moglichen Pfade exponentiell mit der Lange n des HMMs wachst. Zur Losung des Op-
timierungsproblems kommt beispielsweise der sogenannte Viterbi -Algorithmus infrage,
dessen Ziel die Bestimmung des optimalen Pfades (Viterbi-Pfad) bei einer gegebenen
Sequenz der Lange L uber das Alphabet ΣA ist. Viterbi lost das Decodierungsproblem
durch dynamische Programmierung (siehe Kapitel 3.3) mit deutlich weniger Aufwand.
Aufgrund der Einschrankung auf 9 Transitionen pro Modellspalte reduziert sich die
Komplexitat auf O(L · 9n). Fur nahere Details sei auf weiterfuhrende Quellen wie [32],
[21] oder [12] verwiesen.
5.4.3 Alignment mit einem Profil-HMM
Ist das Profil-HMM trainiert und damit die Topologie des HMM festgelegt, so kann
eine Sequenz mit dem HMM aligniert werden. Im Rahmen der Decodierung wird jener
Pfad ermittelt, der die großte Wahrscheinlichkeit oder den hochsten Score aufweist. In
der Abbildung 5.3 wird eine Alignierung der Sequenz d2ez9a2 mit einem aus sieben
Sequenzen bestehenden MSA dargestellt. Modellspalten sind mit dem Buchstaben M
(Match) markiert, Spalten die nicht in das HMM aufgenommen wurden, mit dem
Buchstaben N.
5. Hidden Markov Modelle 61
Abbildung 5.3: Auszug aus einem einfachen Erweiterungs-MSA, bei der den Sequenzendes Profil-HMMs eine einzelne Sequenz hinzugefugt wurde.
Mit den Markierungen m fur Match, d fur Deletion und i fur Insertion werden die
durchlaufenen Zustande des ermittelten Pfads beschrieben. Die Abbildung 5.3 zeigt
den Pfad der mit dem HMM alignierten Sequenz d2ez9a2 als Zeichenkette der Menge
{0,m, i, d,∞}. Wird das Symbole einer Sequenz mit einem m (Match) markiert, so wird
es der nachsten freien Modellspalte M zugeordnet. Wird das Symbol einer Sequenz mit
einem i (Insertion) markiert, so muss das Modell an dieser Position eine zusatzliche
Spalte aufnehmen. Wird das Symbol einer Sequenz mit einem d (Deletion) markiert,
so muss die Sequenz an dieser Position um eine Leerstelle erweitert werden.
Wie Abbildung 5.3 zeigt, kann mit dem aus der Decodierung gewonnenen Pfad eine
Sequenz mit einem MSA aligniert werden. Sollen mehrere Sequenzen aligniert wer-
den, so ist der beschriebenen Prozess mehrmals hintereinander auszufuhren. Wie diese
Aufgabe konkret umgesetzt werden kann, zeigt das folgende Kapitel.
6
Implementierung eines MSA mit
einem Profil-HMM
Das vorliegende Kapitel zeigt die teilweise Umsetzung der in den vorangegangenen
Kapiteln besprochenen Grundlagen und Algorithmen. Die Aufgabe besteht in der Ent-
wicklung, Anwendung und Evaluation eines Verfahrens zur Berechnung von multiplen
Sequenzalignments auf Basis von Profil Hidden Markov Modellen. Als Ausgangsba-
sis der Implementierung wurde das in Zusammenarbeit mit der Universitat Salzburg
(Fachbereich Molekulare Biologie1) und an der Fachhochschule Salzburg (Studiengang
Informationstechnik und System-Management2) entwickelte Softwarepaket HMModeler
[39, 6], eine Werkzeugsammlung zur Erstellung passgenauer Hidden Markov Modelle
fur Proteinfamilien [7, 27, 47], und die damit erstellten Modelle und Daten verwendet.
In einigen Applikationen und Datenstrukturen wurden Anpassungen vorgenommen,
um den geanderten Anforderungen zu genugen. Einige zusatzliche Applikation wurden
ganzlich neu entwickelt. Die Implementierung und das Debugging der Applikationen
und Scripte, sowie die Evaluation der Ergebnisse erfolgte in der im Anhang B.1 doku-
mentierten Umgebung.
Die folgende Kurzbeschreibung der in den nachfolgenden Abschnitten dokumentierten
Implementierungsschritte soll den Ablauf der Entwicklung verdeutlichen.
1Universitat Salzburg, Fachbereich Molekulare Biologie: http://www.uni-salzburg.at/portal/page?_pageid=146,65187&_dad=portal&_schema=PORTAL (Stand: 31. Juli 2010)
2Fachhochschule Salzburg, Studiengang Informationstechnik und System-Management: http://
www.fh-salzburg.ac.at/master/informationstechnologien/(Stand:31.Juli2010)
62
6. Implementierung eines MSA mit einem Profil-HMM 63
6.1 Kurzbeschreibung der Implementierung
Bei der Umsetzung der Problemstellung wird eine mehrstufige und iterative Strategie
verfolgt. Die eigentliche, singulare Alignierung wird mit Hilfe eines Profil-HMM vorge-
nommen, welches auf einem MSA basiert. Mit jeder Alignierung wird eine Reihe von
Scores berechnet. Aus der Menge der Sequenzen werden jene Sequenzen ausgewahlt, die
das hochste Scoring aufweisen und diese dann mit dem Profil-MSA zu einem Ergebnis-
MSA vereint. In einem zweiten Schritt werden aus dem Ergebnis-MSA jene Sequenzbe-
reiche extrahiert, die mit dem Profil-HMM nicht aligniert werden konnen und versucht,
diese Alignierungslucken mit einem alternativen Verfahren zu schließen.
Die folgenden Abschnitte erlautern die schrittweise Umsetzung der Implementierung
und erklaren das Zusammenwirken der zum Einsatz kommenden Methoden und Werk-
zeuge.
6.2 Generieren eines MSA mit einem Profil-HMM
Die Erstellung eines multiplen Sequenzalignments auf Basis eines Profil-HMM wird in
mehrere Teilaufgaben zerlegt. Diese Teilaufgaben werden von einer Reihe von Konso-
lenapplikationen erledigt, welche im Batchbetrieb aufgerufen werden und uber Scripts
miteinander verbunden sind. Der Datenaustausch zwischen den Applikationen erfolgt
in Textdateien, wodurch die Prufung der Teilergebnisse und die Nachvollziehbarkeit
deutlich erleichtert werden.
Folgende Teilaufgaben sind mit den in den Klammern angegebenen Applikationen zu
erledigen:
1. Berechnung der Emissionsfrequenzen auf Basis eines multiplen Sequenzalignments
(recalcParamFile.py)
2. Berechnung der Transitionswahrscheinlichkeiten auf Basis eines multiplen Se-
quenzalignments (transProb.jar)
3. Alignierung aller Sequenzen einer Sequenzdatenbank zum Profil-HMM mit gleich-
zeitigem Scoring der Algnierung (hmmoo.jar)
6. Implementierung eines MSA mit einem Profil-HMM 64
4. Auswahl einer Gruppe von Sequenzen auf Basis der besten Alignment-Scores
(grepseq.jar)
5. Alignierung der ausgewahlten Sequenzen mit dem Profil-MSA zum einem Ver-
einigungs-MSA (amodseq.jar)
Die obige Auflistung der Teilaufgaben entspricht der Reihenfolge der Abarbeitung und
der Reihenfolge der folgenden Beschreibungen.
6.2.1 Training des Profil-HMMs
Die Parametrisierung bzw. das Training des Profil-HMMs erfolgt in zwei Schritten mit
der Berechnung der Emissionsfrequenzen und der Transitionswahrscheinlichkeiten. Da-
tengrundlage des Trainings ist eine Auswahl von Sequenzen mit einer bekannten 3D-
Struktur, die einerseits eine Proteinfamilie umfassend beschreiben, andererseits aber
untereinander keine zu großen Sequenzahnlichkeiten aufweisen sollten. Aus den Se-
quenzen wird automatisch oder manuell ein MSA erstellt. Die Anzahl der Sequenzen
sollte mindestens so groß sein, dass eine reprasentative Schatzung moglich ist (siehe
Abschnitt 5.4.1). Dieses multiple Sequenzalignment bildet die Eingangsdaten zur Be-
rechnung der beiden Wahrscheinlichkeitstypen des Profil-HMM.
6.2.1.1 Berechnung der Emissionswahrscheinlichkeiten
Die Berechnung der Emissionswahrscheinlichkeiten ubernimmt das Python-Script re-
calcParamFile3. Es nimmt ein Parameterfile entgegen und erzeugt unter Angabe zusatz-
licher Parameter, wie beispielsweise der Behandlung der Hintergrundwahrscheinlich-
keiten und eines Gewichtungsfaktors, eine Ausgabedatei. Die Parameterdatei, welche
hier nicht naher beschrieben werden soll, enthalt nebst einigen Verwaltungsdaten im
Wesentlichen eine Beschreibung der alignierten Sequenzen des Profil-MSA und der
Modellspalten. Der Applikationsaufruf erfolgt beispielsweise in der Form:
3Die Applikation recalcParamFile ist Teil des Programmpakets HMModeler und stand fur dieseAufgabenstellung schon zur Verfugung.
6. Implementierung eines MSA mit einem Profil-HMM 65
python recalcParamFile.py -in input.param -out emiss.param ↪→-dirichlet uprior.9compcorr ↪→-alignmentweight 2
Zur Beschreibung der Optionen von recalcParamFile und der eingesetzten Algorithmen
sei auf [6, 39], die Quelltexte und die Applikationhilfe verwiesen. Die Ausgabedatei ent-
spricht in ihrem Grundaufbau weitgehend der Parameterdatei, fuhrt zusatzlich jedoch
eine Matrix mit den Emissionswahrscheinlichkeiten des gegebenen MSA ein. Diese Da-
tei beinhaltet die Eingangsdaten fur den nachsten Berechnungsschritt.
6.2.1.2 Berechnung der Ubergangswahrscheinlichkeiten
Die Berechnung der Ubergangswahrscheinlichkeiten ubernimmt die Java-Applikation
transProbs4. Die Applikation nimmt die vom Script recalcParamFile generierte Er-
gebnisdatei entgegen und erzeugt in Abhangigkeit der Parametrisierung eine HMM-
Datei mit der Beschreibung des Profil-HMM. Nahere Details bzgl. der Berechnungs-
methoden und die Beschreibung der Optionen von transProbs sind in [6] beschrieben.
Der Start der Applikation mit der Einstellung der Pseudocounts zur Transitionsbe-
rechnung (2 2 2) oder der Gewichtung des Expertenwissens (0.1) erfordert folgende
Aufrufsyntax:
java -jar transProbs.jar -p emiss.param -q qas.param ↪→-out current.hmm ↪→-params 2 2 2 0.95 0.7 350 0.1
Die von transProbs erzeugte HMM-Datei enthalt die Beschreibung des Profil-HMMs.
So finden sich darin die Emissions- und Transitionswahrscheinlichkeiten ebenso, wie
die Initial- und Terminalwahrscheinlichkeiten als auch die komplette Beschreibung des
Profil-HMM mit den Signaturen der Modellspalten. Diese Datei dient der Konfiguration
des Modells zur Berechnung erster Alignierungen und Scores.
4Die Applikation transProbs ist Teil des Programmpakets HMModeler und stand fur diese Aufga-benstellung schon zur Verfugung.
6. Implementierung eines MSA mit einem Profil-HMM 66
6.2.2 Single-Scoring und -Alignment mit dem Profil-HMM
Mit der Berechnung der Emissions- und Ubergangswahrscheinlichkeiten ist das Trai-
ning des Profil-HMM abgeschlossen und die Topologie des HMM festgelegt. Die Mo-
dellbeschreibung bildet die statistischen Eigenschaften des Eingangs-MSA soweit ab,
dass damit einzelne Sequenzen zum MSA aligniert und in einem Scoring-Verfahren das
Alignment bewertet werden kann.
Die Alignierung und das Scoring der Sequenzen auf Basis des Profil-HMMs erfolgt
mit der Java-Applikation hmmoo5 [6, 27, 47]. Diese Applikation zeichnet sich dadurch
aus, dass sie einerseits einzelne Sequenzen nach und nach zum Profil-HMM aligniert,
anderseits aber auch eine Reihe von Scores ausgibt, die das Alignment bewerten. In
den nachfolgenden Verfahren werden diese Scores verwendet, um eine Reihung der
Alignierungen durchzufuhren.
Das folgenden Beispiel zeigt den Aufruf der Applikation hmmoo. Die Parameter legen
die HMM-Definitionsdatei und die DB-Datei mit jenen Sequenzen fest, die einzeln
zum Profile-MSA aligniert und mit Scores belegt werden sollen. Der Parameter ’1’
spezifiziert den Datenbanktyp (DBtype Astral6) bzw. das Dateiformat (FASTA [34];
siehe Anhang C.1) und der Parameter ’9’ gibt an, dass alle Scoretypen berechnet werden
sollen (Option Full Scoring). Der Aufruf erfolgt in der Form:
java -jar hmmoo.jar current.hmm 1 astral_scop.fa 9 > hmmoo.results
Der Aufruf der Applikation hmmoo erzeugt eine result-Datei (siehe Anhang C.2), in
der fur jede alignierte Sequenz die Sequenzdaten, die Alignmentinformationen als auch
die berechneten Scores gespeichert sind. Die Alignmentinformationen, wie im folgen-
den Beispiel dargestellt, bestehen aus der Sequenz selbst und aus einer mit der Se-
quenz ausgerichteten Pfadbeschreibung. Diese besteht aus den Zeichen ’M’, ’I’ oder
’D’ und stellt die Zustande dar, die zum Alignment der Sequenz mit dem Profil-MSA
gefuhrt haben. Diese beiden Zeichenketten (siehe folgendes Beispiel) werden in einem
5Die Applikation hmmoo ist Teil des Programmpakets HMModeler und stand fur diese Aufgaben-stellung schon zur Verfugung.
6ASTRAL Compendium for Sequence and Structure Analysis - Databases and Tools: http://
astral.berkeley.edu/ (Stand: 2. Juli 2010)
6. Implementierung eines MSA mit einem Profil-HMM 67
spateren Arbeitsschritt (Abschnitt 6.2.4 diese Kapitels) zur Berechnung eines multi-
plen Alignments ausgewertet. Eine Beispielsequenz mit Pfadbeschreibung in Form einer
MID-Zustandsmaske lautet:
MSA:
SLFEQLGGQAAVQAVTAQFYANIQADATVATFFNGIDMPNQTNKTAAFLCAALGGPNAWTSLFEQLG ↪→GQAAVQAVTAQFYANIQADATVATFFNGIDMPNQTNKTAAFLCAALGGPNAWTGRNLKEVHANMGVS ↪→NAQFTTVIGHLRSALTGAGVAAALVEQTVAVAETVRGDVVTV
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIMM ↪→MMMMMMMMMMMMMMMMMMMMMMMMMMMIIMMMMMMMMMMMMMMMMMMMIMMIMMMMMMMMMIMMMMM ↪→MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMIIIIIIIIII
Wird die Option ’Full Scoring ’ gewahlt, so werden zwolf unterschiedliche Scores berech-
net, welche die Qualitat des Alignments zwischen dem Profil-MSA und der jeweiligen
alignierten Sequenz bewerten. Ein Scoring-Block, wie im folgenden Auszug eines result-
Datensatzes fur die Sequenz d1dlwa aus der Familie a.1.1.1 dargestellt, besteht nebst
oben beschriebenen Informationen aus der jeweiligen Score-Bezeichnung und mit ei-
nem ’:’ davon getrennt (Delimiter) aus dem Score selbst, der als Fließkommazahl im
Double-Format ausgewiesen wird.
Plain Score: -576.7240977002209
Plain Score (Forward): -569.0899435594409
Simple Score: -647.0919876918512
Simple Score (Forward): -647.0919876918512
Reversed Score: -581.9745888906357
Reversed Score (Forward): -572.879493576253
Simple Corrected Score: 70.36788999163025
Simple Corrected Score (Forward): 78.00204413241022
Reverse Corrected Score: 5.2504911904147775
Reverse Corrected Score (Forward): 3.789550016812086
HMMer Corrected Score: 65.11109477813987
HMMer Corrected Score (Forward): 71.77954943225613
Welche Scores von hmmoo berechnet werden, wie diese zu interpretieren sind und auf
Basis welcher mathematischen Grundlagen sie wie berechnet werden, ist Teil der Doku-
mentation in [6]. Festzuhalten ist, dass die Großen der Scores einzeln betrachtet kaum
aussagekraftig sind. Erst der Vergleich der Scores untereinander und eine Reihung
der Scores lasst eine Aussage daruber zu, ob das Scoring einer Sequenz ein”gutes“,
6. Implementierung eines MSA mit einem Profil-HMM 68
”maßiges“ oder gar
”schlechtes“ Alignment anzeigt. Eine Auswahl der Sequenzen auf-
grund der Reihung der Scores ubernimmt die im Folgenden beschriebene Applikation
grepseq.
6.2.3 Auswahl der zu alignierenden Sequenzen
Die fur diese Arbeit entwickelte Java-Applikation grepseq liest die von hmmoo generier-
te result-Datei und ordnet die Reihen nach vorgegebenen Kriterien. Die Applikation
grepseq bietet die Option7, bestimmte Datensatze uber regulare Ausdrucke (regular
expressions) explizit einer Verarbeitung zuzufuhren oder sie auch ausdrucklich auszu-
schließen. Ebenso konnen die Datensatze auf ihre Scores hin untersucht und nach den
Scores gereiht werden, bevor aus den verbliebenen Datensatzen jene selektiert werden,
welche die hochsten Scores erhalten haben.
Die einfachste Moglichkeit der Sequenzauswahl ist eine absteigende Sortierung der Sco-
res, um aus allen Datensatzen jene n Datensatze auszuwahlen, welche bei einem Sco-
retyp die hochste Bewertung erhalten haben. Diese Methode ware legitim, zeigt aber
in der Praxis einige wesentliche Nachteile. Kommt es beispielsweise zu Ausreißern in
der Bewertung der Alignierung, das heißt zu falschen und besonders niedrigen oder
hohen Scores, so wurden genau diese Scores in einem spateren Schritt das Bewertungs-
ergebnis verfalschen. Wird das MSA zu einem spateren Zeitpunkt zum Training eines
Profil-HMM weiterverwendet, so sollen nur jene Sequenzen in das MSA aufgenommen
werden, die den Sequenzen im Profil-MSA nicht zu ahnlich sind. Gleichzeitig darf als
limitierende Schranke aber die Spezifitat des HMMs nicht zu gering werden. Es durfen
demnach nur Sequenzen mit hinreichender Ahnlichkeit in das HMM aufgenommen
werden. Grundlage dieser Uberlegung ist die Feststellung, dass aus informationstheo-
retischer Sicht der Informationsgehalt eines MSA sinkt, wenn viele nahezu gleiche Se-
quenzen das MSA beschreiben. Dagegen steigt der Informationsgehalt des MSA, wenn
das MSA aus moglichst verschiedenen Sequenzen gebildet wird, obwohl diese eine ho-
he strukturelle Ahnlichkeit aufweisen. Gemessen wird der Informationsgehalt mit der
7Eine Kurzbeschreibung der Optionen wird ausgegeben, wenn die Applikation mit der Option -hoder –help aufgerufen wird. Siehe Anhang B.2.
6. Implementierung eines MSA mit einem Profil-HMM 69
Entropie, die zu einem spateren Zeitpunkt in dieser Arbeit noch einmal diskutiert und
berechnet wird.
Die Losung beider Probleme besteht darin, Ausreißer und Sequenzen mit auffalligen
Scores auszuschließen und nur jene weiterzuverarbeiten, die obige Kriterien erfullen.
Diese Losung bedarf der genauen Festlegung des Begriffs Ausreißer und der Klarung
der Frage, wie eine Variable als Ausreißer erkannt werden kann, wenn der Wertebereich
nicht bekannt ist. Dafur wurden in grepseq eine Losung implementiert, der statistische
Uberlegungen zugrunde liegen:
Untersucht man die Scores nach deren Verteilung, so stellt man fest, dass diese in
den meisten Fallen einer Verteilungsfunktion entsprechen. So zeigt beispielsweise die
Untersuchung des”Reverse Corrected Score“ von etwa 9500 Datensatzen, dass dieser
annahernd einer Normalverteilung entspricht (siehe Abbildung 6.1). Dieser Erkenntnis
entsprechend kann aus dem Mittelwert µ und der Standardabweichung σ der Scores
eine Normalverteilung angenahert werden.
Abbildung 6.1: Reale Verteilung und Annaherung der Normalverteilung fur ReverseCorrected Scores auf Basis des Mittelwerts fur Sequenzen der Familie d.38.1.5.
Mit der Annahme, ein Score wurde zufallig sein und die Verteilung eines Scores in
allen Datensatzen wurde einer Normalverteilung entsprechen, kann uber die kumulati-
ve Verteilungsfunktion (cumulative distribution function; CDF ) die Wahrscheinlichkeit
berechnet werden, mit der ein Score einen Wert in einem bestimmten Bereich einnimmt.
Realisiert wurde dies uber die Moglichkeit der Angabe eines Grenzwertes (Threshold)
6. Implementierung eines MSA mit einem Profil-HMM 70
0 <= t <= 1 , der einen Verwertungsbereich markiert (siehe Abbildung 6.1; rot mar-
kierter Bereich rechts des Thresholds). Liegt ein Score im Ablehnbereich, also links des
Thresholds (s < t), so wird der Datensatz nicht mehr weiter beachtet. Liegt der Score
auf dem oder rechts des Thresholds (s >= t), also im Verwertungsbereich, so wird die
Sequenz weiterverarbeitet.
Problematisch ist bei der Schatzung des Mittelwerts, dass extreme Ausreißer den Mit-
telwert deutlich verfalschen konnen. Diese Uberlegung wurde in grepseq dahingehend
verwertet, dass von der Berechnung des Mittelwerts alle jene Sequenzen ausgeschlossen
werden konnen, deren Wert derart im Randbereich der Verteilung liegt, dass es nicht
mehr glaubhaft ist, ein derartiges Ergebnis zu erhalten. Ebenso lasst sich dieser Um-
stand damit umgehen, dass anstatt des Mittelwerts der Median berechnet wird und
die Normalverteilung um den Median positioniert wird. Gerade bei Normalverteilungen
liegt nahe, dass im Bereich des Mittelwerts ein großer Teil der Scores liegen sollte. Der
Abdrift des Medians durch extreme Ausreißer ist geringer als beim Mittelwert. Dieser
Umstand wurde in grepseq entsprechend berucksichtigt.
Eine bei den Testlaufen haufig verwendete und im folgenden Beispiel dargestellte Auf-
rufsequenz der Applikation grepseq legt fest, dass die Sortierung nach dem Reverse
Corrected Score (RCS) erfolgt. Weiters wird die NormCDF des RCS uber den Median
(NCDFM) berechnet und nur jene Sequenzen selektiert, deren Scores uber oder gleich
(greater equal - GE) dem angegebenen Threshold von 0.999 liegen.
java -jar grepseq.jar -s RCS ↪→-c RCS NCDFM -t GE 0.999 -n 3 ↪→-o hmmoo.results topseqs.results
Die Applikation grepseq liest alle Sequenzen E := {SEQ1, SEQ2, · · ·SEQu} aus der
Datei hmmoo.results und wahlt jene Sequenzen aus, deren Score jenseits des frei
wahlbaren Threshold t liegen. Wie sich in der Bewertung der Alignments (siehe Ka-
pitel 7) zeigt, gilt ein Threshold von 0.999 als großzugig gewahlt. Aus diesem Grund
fallen unter bestimmten Umstanden zu viele Sequenzen in den durch den Threshold
begrenzten Verwertungsbereich. Liegen zu viele Sequenzen oberhalb der angegebenen
Grenze, so kann die Anzahl der Sequenzen begrenzt werden. Im Beispiel wurde die
Auswahl auf jene Sequenzen mit den 3 hochsten Scores limitiert.
6. Implementierung eines MSA mit einem Profil-HMM 71
Nachdem die Sequenzen der Eingangsmenge E := {SEQ1, SEQ2, · · ·SEQu} einzeln
mit dem Profil-MSA MSAP := {SEQ1, SEQ2, · · · SEQm} aligniert und mit grepseq
gefiltert wurden, sind nur noch jene Sequenzen ubriggeblieben, die den Erfolg verspre-
chendsten Score aufweisen. Sie bilden eine stark verkleinerte Menge V aus den mit dem
Profil-MSA voralignierten Sequenzen.
6.2.4 Progressives Alignment mit dem Profil-HMM
Die Aufgabe des folgenden Prozessschritts ist die Vereinigung der voralignierten Se-
quenzen V mit dem Profil-MSA MSAP zu einem neuen multiplen Sequenzalignment
mit den Sequenzen MSAV = MSAP ∪ V .
Die implementierte Methode ist mit jener der im Kapitel 4.3 eingefuhrten iterativen
und progressiven Vorgangsweise vergleichbar, unterscheidet sich jedoch darin, dass die
Reihenfolge der Alignierung nicht uber einen binaren Baum bestimmt wird. Stattdessen
werden in jedem Iterationsschritt j = {1, 2, · · · |V |} die Sequenzen SEQi in der von
grepseq festgelegten Reihenfolge schrittweise zum MSA aligniert. Grepseq ordnet die
Sequenzen dazu in der Reihenfolge ihrer Scores s fur das globale Alignment von SEQi
mitMSAP , folglich werden die Sequenzen mit den hochsten Scores am Beginn, und jene
mit niedrigsten Scores am Ende des Vereinigungsprozesse dem Profil-MSA hinzugefugt.
Die Aufgabe der multiplen Alignierung ubernimmt die im Rahmen dieser Arbeit ent-
wickelte Java-Applikation amodseq. Amodseq liest die Profil-HMM-Parameterdatei mit
dem Profil-MSA MSAP und die von grepseq erzeugte Datei mit den vorausgewahlten
Sequenzen V . Im ersten Schritt j = 1 wird jene Sequenz aus V mit dem hochsten Sco-
re dem bestehenden Profil-MSA MSAP hinzugefugt und ein neues Vereinigungs-MSA
(VMSA) MSAj gebildet. Im ersten Schritt entsteht das MSA1.
x1 = argmaxi
si
MSA1 := MSAP ∪ {SEQx1}(6.1)
6. Implementierung eines MSA mit einem Profil-HMM 72
In jedem folgenden Iterationsschritt j → j+ 1 wird eine weitere Sequenz SEQj ∈ V in
absteigender Reihenfolge der Scores dem bestehenden MSAj−1 hinzugefugt und daraus
wiederum ein neues MSAj berechnet.
xj = argmaxi 6=x1,...,xj−1
sj
MSAj := MSAj−1 ∪ {SEQxj}
(6.2)
Dieser Vorgang muss derart implementiert werden, dass sowohl die Alignierungen des
aktuellen MSAj−1, als auch die der aktuell hinzuzufugenden Sequenz SEQxjmit dem
ursprunglichen Profil-MSA MSAP gultig bleiben. Folgende Situationen mussen dafur
beim Alignment einer Sequenz SEQxjmit den multiplen Sequenzalignment MSAj−1
berucksichtigt werden:
• Falls die mit dem MSAj−1 alignierte Sequenz SEQxjan der Position k keine
weitere Lucke in MSAj−1 einfugt (das heißt den Zustand ’D’ oder ’M’ signalisiert),
so kann das Symbol an der Position k in das MSAj ubernommen werden.
• Falls die mit dem MSAj−1 alignierte Sequenz SEQxjan der Position k eine
zusatzliche Lucke einfugt, so muss diese Lucke in samtliche Sequenzen der MSAj
eingefugt werden, um die Sequenz dem MSAj−1 hinzufugen zu konnen.
Die korrekte Behandlung dieser Situationen stellt eine der Kernfunktionen der Appli-
kation amodseq dar. Jene Codeteile, welche die Alignierung vornehmen und dabei die
oben beschriebenen Anforderungen berucksichtigen, werden im Anhang D, Listing D.1
auszugsweise gezeigt.
Am Ende dieses Prozessschritts erzeugt amodseq eine Datei mit der Vereinigungs-MSA.
Eine Prufung der Ergebnisse ergibt, dass die Alignierung der Sequenzen mit dem Profil-
MSA funktioniert (eine detaillierte Prufung der Vereinigungs-MSA erfolgt in Kapitel 7),
jedoch augenscheinlich einige Schwachen zeigt, die im folgenden Abschnitt Beachtung
finden.
6. Implementierung eines MSA mit einem Profil-HMM 73
6.2.5 Kombination mehrerer Alignmentmethoden
Die Sichtkontrolle des von amodgrep generierten Vereinigungs-MSA ergibt, dass die
Alignierung mit dem Profil-HMM plausible Ergebnisse ausweist. Eine genauere Prufung
der Resultate zeigt jedoch auch, dass unter bestimmten Voraussetzungen bestimmte
Sequenzregionen ganzlich unaligniert bleiben. Die Abbildung 6.2 zeigt ein Profil-MSA
bestehend aus vier Sequenzen, zu dem drei weitere Sequenzen aligniert wurden. Die vier
alignierten Profil-Sequenzen (in der Abbildung werden nur Ausschnitte von Position
20 bis 60 dargestellt), spezifizieren das Profil-HMM und legen in der Trainingsphase
die Topologie des HMMs fest.
Abbildung 6.2: Auszug aus einem unvollstandigen MSA auf Basis eines Profil-HMM.
Jene Spalten des Profil-MSA, die mit dem Zustand M (Match) markiert sind, sind
Teil des Profil-HMM, da sie im Training des HMM einen Konsens aufweisen. Diese
Regionen wurden vom Profil-HMM zur Alignierung herangezogen (grune Bereiche).
Jene Sequenzbereiche, die im Profil-MSA keinen Konsens aufweisen (weiße Spalten
mit dem Zustand N), wurden nicht in das Profil-HMM aufgenommen und bleiben
bei der Alignierung weitgehend unberucksichtigt. Die Folge davon ist, dass bestimmte
Regionen (fett umrandete orange Bereiche) der hinzugefugten Sequenzen unaligniert
bleiben.
Die unalignierten Teile der Sequenzen werden bei jedem Interationsschritt in das ak-
tuelle MSA MSAi eingefugt, was zur Folge hat, dass sich alle Sequenzen des aktuellen
MSAs entsprechend verlangern. Sichtbar wird dieser Effekt im Beispiel der Abbildung
6.2 in den Spalten 49 bis 55, in denen eine Lucke im MSA entstanden ist, die innerhalb
6. Implementierung eines MSA mit einem Profil-HMM 74
von drei Interationsschritten angewachsen ist. Dabei wird deutlich, welche Auswir-
kungen die Alignierung mehrerer Sequenzen mit deutlich langeren Lucken zwischen
den Match-Bereichen auf das Vereinigungs-MSA hat. Das VMSA wachst in der Breite
uberproportional an.
Eine mogliche Gegenmaßnahme besteht darin, diese entstandenen Lucken uber alter-
native Verfahren zu alignieren. Unter Anwendung von PAM- oder BLOSUM-Matrizen
ist eine nachtragliche Alignierung dieser Bereiche durchaus denkbar. Die Applikation
amodseq bedient sich jedoch des ClustalW-Algorithmus um diese Lucken zu schließen.
Zu diesem Zweck werden nach der Erstellung des Vereinigungs-MSA jene Bereiche
aus dem MSA extrahiert (in der Abbildung 6.2 sind dies die umrandeten weißen und
orangefarbigen Bereiche der Spalten N und i), die vom HMM nicht aligniert werden.
Es entstehen sogenannte MSA-Snippets oder Profile8. Dabei wird zwischen den Profil-
(Profil1 ; weiße Bereiche) und den hinzugefugten Bereichen (Profil2 ; orange Bereiche)
unterschieden. Zu diesem Zeitpunkt kommt die Applikation ClustalW29 zum Einsatz,
um die beiden Profile zu alignieren, wobei die Alignierung des Profil1 unverandert
bleibt und das Profil2 passend dazu aligniert wird. Dieser Arbeitsschritt wird fur alle
Lucken bzw. Snippets wiederholt. Am Ende werden die von ClustalW alignierten MSA-
Snippets wieder in das Vereinigungs-MSA zuruckgeschrieben. Die folgende Abbildung
6.3 zeigt das Ergebnis dieses Ablaufs.
Deutlich sichtbar wird der Erfolg dieser Methode, vergleicht man die Breite der Lucken
nach der ClustalW-Alignierung mit den Lucken zuvor. Der Bereich der Spalten 43-45
(vormals die Spalten 49 bis 55) kann mit ClustalW beispielsweise erfolgreich aligniert
und die Lucke damit deutlich verkleinert werden. Das Vereinigungs-MSA wird dadurch
kompakter bzw. dessen Sequenzen entsprechend verkurzt.
Der Aufruf der Java-Applikation amodseq fur das oben beschriebene Szenario lautet:
java -jar amodseq.jar -p current.hmm -a topseqs.results ↪→-o current.msa -d ↪→-x 3 -TYPE=PROTEIN -GAPOPEN=10 -GAPEXT=0.1
8Der Begriff Profile wird alternativ verwendet, um damit der Terminologie von ClustalW zu folgen.ClustalW werden diese Snippets als Profile1 und Profile2 zugefuhrt, um diese zu alignieren.
9CLUSTAL 2.0.12 Multiple Sequence Alignments; http://www.clustal.org/ (Stand: 13. April2010)
6. Implementierung eines MSA mit einem Profil-HMM 75
Abbildung 6.3: Auszug aus einem MSA auf Basis eines Profil-HMM und ClustalW
Die Parameter -x 3 legen fest, dass die darauf folgenden drei Parameter unverandert
an ClustalW2 durchgeschliffen werden. Im obigen Beispiel werden die Gap-opening-
und Gap-extension Penalties mit den Werten 10 und 0.1 belegt. Eine Reihe weite-
re Optionen10 bieten zusatzliche Optimierungsmoglichkeiten im Zusammenwirken von
amodseq mit der Applikation ClustalW11. Einige Parameter, die im Zusammenhang
mit ClustalW immer zur Anwendung kommen, sind in der Konfigurationsdatei amod-
seq.ini abgelegt (siehe Anhang B.4).
6.2.6 Progressive Expansion eines Profil-HMM
Wie aus der Einfuhrung im Kapitel 4 hervorgeht, dient die Erstellung eines multiplen
Sequenzalignments in erster Linie der Abbildung evolutionarer Zusammenhange in der
Entwicklung von Sequenzen. Gleichzeitig hat der Abschnitt 5.4.1 gezeigt, wie sich aus
multiplen Sequenzalignments Profile erstellen lassen, deren Aufgabe es ist, ein MSA so
zu beschreiben, dass damit strukturell verwandte Sequenzen gefunden bzw. aligniert
werden konnen.
Ausgangspunkt des bisher beschriebenen profilbasierten Verfahrens ist ein manuell oder
ein mit Computerapplikationen automatisch erstelltes multiples Alignment aus oft nur
10Eine Kurzbeschreibung der Optionen von amodseq wird ausgegeben, wenn die Applikation mitder Option -h oder –help aufgerufen wird. Siehe Anhang B.3.
11
”Using ClustalX for multiple sequence alignment“ ist unter http://www.clustal.org/download/
clustalw_help.txt verfugbar. Ein Aufruf der Applikation Clustal2 mit dem Parameter -HELP zeigtebenfalls die von ClustalW zur Verfugung stehenden Optionen.
6. Implementierung eines MSA mit einem Profil-HMM 76
einigen wenigen Sequenzen. Dieses Profil-MSA bestimmt im Training die Topologie
des HMMs. In der Regel wird ein MSA vorbereitet, welches eine bestimmte Fami-
lie moglichst umfassend und exakt beschreibt. Die Grundlage dieser Beschreibungen
bilden im Allgemeinen jene Sequenzen, deren strukturelle Ahnlichkeiten und deren
evolutionare Verwandtschaft hinlanglich bekannt und uberpruft sind.
Erstellt man aus einem Profil-MSA und einigen Sequenzen ein neues MSA, so liegt
die Idee nahe, die expandierten MSAs als Profil-MSAs wiederzuverwenden und uber
mehrere Zyklen hinweg damit wiederholt einen neuen Alignierungsprozess zu starten.
Dieser progressive Ansatz fuhrt dazu, dass sich die HMM-Topologie bei jedem Zyklus
dynamisch verandert und an das neue Profil anpasst.
Die symbolische Darstellung in Abbildung 6.4 beschreibt das Prinzip dieser Methode.
Rx und Ry symbolisieren die Achsen eines gedachten Raums zur Darstellung der Scores,
wobei damit alle moglichen durch das HMM aus einer Sequenz ableitbaren Informatio-
nen verstanden werden. Die Sequenzen werden in diesem Informationsraum als Punkte
eingezeichnet. Der Informationsgehalt des Profil-HMM wird in der Abbildung 6.4 als
die vom Score-Threshold eingegrenzte Flache symbolisiert. Die Modellmittelpunkte
sind als Schwerpunkt der Score-Threshold-Flache oder Sequenzpunkte vorstellbar und
mit einem Kreuz eingezeichnet.
Abbildung 6.4: Symbolische Darstellung der progressiven HMM-Expansion.
6. Implementierung eines MSA mit einem Profil-HMM 77
Die in der Abbildung 6.4 beschriebene Ausgangssituation entspricht dem (hellblauen)
Startmodell, das mit neun (roten) Sequenzen beschrieben wird. Im ersten Alignierungs-
zyklus werden dem Profil-MSA funf (grune) neue Sequenzen hinzugefugt und daraus
fur einen weiteren Zyklus ein neues (hellrotes) Profil-MSA mit 14 Sequenzen erzeugt.
Im zweiten Zyklus werden mit dem Profil-MSA vier weitere (blaue) Sequenzen gefun-
den, die aufgrund der Expansion des HMMs nun in die Nahe des Profil-MSA geruckt
sind. Der eingezeichnete Pfad der Modellmittelpunkte verdeutlicht den Drift des Mo-
dells im Laufe des Prozesses. Das expandierte Profil-HMM andert durch die Aufnahme
neuer Sequenzen nicht nur seinen Informationsgehalt, sondern auch seinen Modellmit-
telpunkt. Expandiert in jedem Zyklus das Modell, so konnen im Idealfall damit in den
Folgezyklen Sequenzen gefunden werden, die mit dem Ausgangsmodell (noch) nicht
gefunden wurden.
Die praktische Umsetzung dieses Verfahrens erfordert Werkzeuge, um aus der ent-
standenen Alignierung neue Parameterdateien zu erstellen. Die dafur entwickelte Java-
Applikation modmsapm.jar erfullt den geforderten Zweck und liest die aktuelle param-
Datei, mit dem die letzte Alignierung durchgefuhrt wurde, und die von amodseq ge-
nerierte MSA-Beschreibung, um daraus eine neue Parameterdatei fur den nachsten
Zyklus zu erstellen:
java -jar modmsapm.jar -p input.param -m current.msa ↪→-o nextinput.param
Wesentlich bei der aktuellen Implementierung ist die Einschrankung, dass keine neu-
en Modellspalten hinzugenommen werden, sondern das Profil-MSA ausschließlich um
zusatzliche Sequenzen erweitert wird. Das Einfugen zusatzlicher Modellspalten ist mit
wenig technischen Aufwand verbunden, hatte aufgrund der komplexen Auswirkungen
den Vergleich der Ergebnisse aber deutlich erschwert und den Aufwand der Evaluation
(siehe nachster Abschnitt) deutlich erhoht.
Getaktet werden die Zyklen uber ein Bash-Script, welches eine feste Anzahl von Zyklen
implementiert. In jedem dieser Zyklen werden die oben genannten Dateien erzeugt und
alle Dateien zum Zweck einer spateren Evaluation gesichert. Zusatzlich werden einige
6. Implementierung eines MSA mit einem Profil-HMM 78
Dateien mit Zwischenergebnissen generiert, die eine spatere Auswertung unter Zuhilfe-
nahme zusatzlicher Scripts und externen Werkzeugen (zum Beispiel mit MATLAB12)
deutlich erleichtern.
Alle Scriptdateien, Auszuge aus den Quelltexten der Applikationen, die in dieser Arbeit
verwendeten Testdaten und daraus generierte Ergebnisdateien sind auf dem beigefugten
Datentrager gesichert (siehe Anhang E).
12MATLAB ist eine Produktfamilie der Firma The MathWorks, http://www.mathworks.com/
products/matlab/ (Stand: 24. August 2010)
7
Bewertung der Ergebnisse
Die bisherigen Ausfuhrungen haben die Implementierung beschrieben, die objektive
Prufung der Ergebnisse ist bislang ausgeblieben. Dieser Abschnitt beschreibt aus-
gewahlte Methoden zur Evaluation der Ergebnisse und pruft, inwieweit die Metho-
de der progressiven Expansion des Profil-HMM zur schrittweisen Alignierung genutzt
werden kann.
7.1 Testdaten und -settings
Die Bewertung der Ergebnisse der Verfahren und Applikationen erfolgt nach Testlaufen
unter folgenden Rahmenbedingungen:
• Die Sequenzdatenbank astral_scop_40.1.73 enthalt 9536 Sequenzen aus un-
terschiedlichsten Familien und mit unterschiedlichen Sequenzlangen.
• Das Profil-HMM wurde mit 20 MSAs ausgewahlter Familien trainiert. Die
Parameterdateien enthalten neben den Profil-MSAs auch zusatzliche Angaben
daruber, welche Spalten als Modellspalten ausgewahlt werden; dafur ist auch
Expertenwissen eingeflossen, um moglichst aussagekraftige und korrekte Profil-
MSA zu erhalten.
• Ein kompletter Testlauf bestand aus 4 bis 6 Zyklen.
• Die Aufnahme einer neuen Sequenz wurde uber den CDF-Threshold von 0.999
begrenzt; nach oben hin fand keine Begrenzung statt.
79
7. Bewertung der Ergebnisse 80
• Das Profil-MSA wurde in jedem Zyklus um maximal 3 Sequenzen erweitert.
Jedes der in der Tabelle 7.1 angefuhrten multipler Sequenzalignments zum Training
eines Profil-HMMs wurde aus den alignierten Sequenzen einer bestimmten Familie ge-
bildet. In der Tabelle sind sowohl die Familiennamen der Profil-MSA-Sequenzen als
auch die Anzahl der Spalten der Profil-MSA und die Anzahl der im Profil-HMM ent-
haltenen Emissionsspalten angegeben.
Tabelle 7.1: Multiple Sequenzalignments zur Beschreibung von Profil-HMMs. Die in derSpalte ’MSACols’ angegebenen Werte beschreiben die Lange der Sequenzen im Profil-MSA, die Spalte ’Col.Emiss’ beschreibt die Anzahl der Emissionsspalten im HMM.Die Werte der Spalten ’Modell-Sequenzen’ von ’Z1’ bis ’Z6’ nennen die Anzahl derSequenzen, aus denen im jeweiligen Zyklus das Profil-MSA gebildet wurde. Der Test indiesem Beispiel fuhrte uber sechs Zyklen, in dem in jedem Zyklus drei neue Sequenzenzum Profil-MSA hinzugefugt wurden.
7. Bewertung der Ergebnisse 81
Im rechten Teil der Tabelle ist angefuhrt, um wie viele Sequenzen das Modell pro
Zyklus erweitert wurde. Wie die Tabelle zeigt, wurde bei allen Familien bei jedem
Durchlauf drei Sequenzen gefunden, die uber dem vorgegebenen CDF-Threshold lagen.
Damit wurde bei jedem Zyklus die Profil-MSAs um die maximale Anzahl der erlaubten
Sequenzen erweitert.
Jede der 9536 Sequenzen aus der im FASTA-Format vorliegenden Testdatenbank ist
einer Familie zugeordnet. Die Sequenzen haben unterschiedliche Langen im Bereich
von 21 bis 1504 Zeichen. Die Verteilung der Sequenzlangen ist unregelmaßig und un-
abhangig von deren Familienzugehorigkeit. Wie die Sequenzlangenverteilung in Abbil-
dung 7.1 zeigt, hat der Hauptanteil der Sequenzen eine Lange im Bereich von 50 bis
180 Zeichen.
Abbildung 7.1: Langenverteilung aller Sequenzen aus der Testdatenbank’astral scop 40.1.73’
Wie im Kapitel 6.2.3 beschrieben, aligniert die Applikation hmmoo nicht nur ein-
zelne Sequenzen zum Profil-MSA, sondern bewertet gleichzeitig die Alignierung mit
einer Reihe von Scores. Wie in der Methode der Vorauswahl der Sequenzen schon
berucksichtigt wurde, folgen die Scores bestimmten Verteilungen, die fur die Zwecke
dieser Arbeit mit einer Normalverteilung angenahert werden. Aus diesem Grund bil-
den sowohl die Langenabhangigkeiten der Scores als auch die Dichtefunktionen (Wahr-
scheinlichkeitsverteilung) der Scores die Grundlage der Evaluation der multiplen Se-
quenzalignments, wie im folgen Abschnitt dargelegt wird.
7. Bewertung der Ergebnisse 82
7.2 Dichtefunktion der Scores
Auer et al. [6, 27] haben die Applikation hmmoo und die Verteilung und Eigenschaften
der Scores hinreichend dokumentiert. Sie haben in ihrer Arbeit aber ebenso festgehal-
ten, dass nicht alle Scores gleichermaßen geeignet sind, die Alignierung einer Sequenz
mit einem Profil-MSA uber Scores zu beurteilen. So unterstutzt die Applikation hmmoo
wohl 12 verschiedene Scoretypen, einige davon sind aber stark abhangig von der Anzahl
der Emissionsspalten und/oder von der Lange der Sequenzen. Damit ist ein Vergleich
gleicher Scoretypen uber Sequenzen - und weniger noch uber Profil-MSAs hinweg -
nur schwer oder nicht moglich. Durch die Einfuhrung von Corrected Scores konnte die
Langenabhangigkeit zwar deutlich reduziert, jedoch nicht ganzlich korrigiert werden.
Aligniert man alle Sequenzen der Sequenzdatenbank astral_scop_40.1.73 mit ei-
nem Profil-HMM, so erhalt man 9536 Scores pro Scoretyp. Wird aus den Scores eines
Typs die Wahrscheinlichkeitsdichtefunktion (propability density function) errechnet, so
erhalt man ein Bild der wahrscheinlichen Verteilung der Scores uber den Wertebe-
reich des jeweiligen Scoretyps. Diese Darstellungsform lasst sich auch dafur nutzen, die
Trennscharfe eines Modells fur die jeweiligen Scores sichtbar zu machen. Dazu wer-
den die Scores nach Family, Superfamily, Fold und Other (Other = Class + Rest der
Sequenzen; siehe Abschnitt 2.5.1) getrennt ausgewertet und daraus die Wahrschein-
lichkeitsdichtefunktion errechnet [47].
Abbildung 7.2 zeigt die Dichtefunktionen fur den Reverse Corrected Score1 (RCS) der
Profil-MSA Familie b.121.4.1. Wie aus dem Diagramm deutlich wird, setzen sich die
Scores jener Sequenzen, die aus der gleichen Familie wie jene des Profil-MSA stammen,
deutlich vom Rest der Scores ab. Sequenzen der Familie b.121.4.1 werden im Mittel
etwa mit dem Score 68 bewertet. Erhalt man einen Score von 60, so handelt es sich
mit hoher Wahrscheinlichkeit um eine Sequenz der gleichen Familie, zumal der Rest
der Sequenzen deutlich unter 60 bewertet wurde. Selbst der Score von 40 ist noch sehr
eindeutig.
Die restlichen Sequenzen unterscheiden sich bezuglich deren RCS kaum voneinander.
Sowohl die Sequenzen der Fold b.121 als auch die Sequenzen ganz anderer Klassen
1Zur detaillierten Erklarung der Scoretypen sei auf die Arbeit von Auer [6] verwiesen.
7. Bewertung der Ergebnisse 83
Abbildung 7.2: Dichtefunktionen des Reverse Corrected Scores fur das Profil-MSA derFamilie b.121.4.1
werden im Mittel mit dem Score um Null bewertet. Einzig die Sequenzen der gleichen
Superfamily wie das Profil-HMM, werden mit einem hoheren Scoring bewertet und
heben sich, wenn auch nur minimal, ab. Die Dichtekurve zeigt aber auch, dass der
Abstand der Scores der Superfamily-Sequenzen zu den restlichen Sequenzen zu gering
ist, um Superfamilies von den Folds und Others sicher zu unterscheiden.
Wie sollte eine Dichtefunktion optimalerweise aussehen? Hierbei entscheidet weni-
ger die Form der Funktion, also vielmehr die Lage der Dichtefunktionen zueinander.
Gewunscht ist eine Verteilung der Scores, so dass aufgrund des Scores schon eine mit
hoher Wahrscheinlichkeit korrekte Klassifizierung erfolgen kann. Die Scorebereiche, im
Diagramm als die Breite der Dichtefunktionen sichtbar, sollten sich demnach moglichst
wenig uberlappen.
Berechnet man die Scores einzelner Sequenzen verglichen uber die verschiedenen Mo-
delle, so fallt auf, dass einige der Scores zwar deutlich weniger abhangig von der Lange
der Sequenzen sind als andere2 [6], trotzdem aber noch Abhangigkeiten zwischen dem
Modell und den Sequenzen bestehen, die den Score entsprechend anheben oder senken
konnen. Berechnungen zeigen, dass unter bestimmten Voraussetzungen die Varianz der
Scores groß wird, in anderen Situation aber wiederum vergleichsweise klein. Versucht
man nun, die Bewertung der Klassifizierungsfahigkeit eines Modells vom Abstand des
Mittelwerts der Dichtefunktionen abhangig zu machen, so kann keinesfalls auf absolute
2Diese langenunabhangigen Scores werden von der Applikation hmmo immer mit dem Zusatz Cor-rected versehen.
7. Bewertung der Ergebnisse 84
Großen zugegriffen werden. Dichtefunktionen die eine große Standardabweichung anzei-
gen, mussen einen deutlich großeren Abstand untereinander einnehmen, um sich nicht
zu uberlappen, als Dichtefunktionen mit kleineren Standardabweichungen. Die Große
die es zu bewerten gilt ist vielmehr der relative Abstand der Dichtefunktionen unter-
einander. Diese Uberlegung macht eine vorherige Standardisierung der Scores sinnvoll.
7.3 Standardisierung der Dichtefunktion
Um die Wahrscheinlichkeitsverteilungen untereinander vergleichbar zu machen, muss
eine allgemeine Verteilung in eine standardisierte Verteilung transformiert werden. Eine
allgemeine Gaußsche Normalverteilung S ∼ N(µ, σ2) mit den Parametern µ und σ lasst
sich in eine Standardnormalverteilung uberfuhren, in der der Mittelwert µ = 0 und die
Standardabweichung σ = 1 wird. Hierfur werden normalverteilte Scores S durch eine
lineare Transformation zu sog. standardisierten Scores S ∼ N(0, 1). Die Umrechnung
erfolgt gleich der Formel der Z-Scores (siehe Gleichung 3.10) in der Form [37]:
S =S − µσ
(7.1)
Wird diese Standardisierung auf die Scores aller Sequenzen angewendet, indem die
Scores mit dem Mittelwert µO und der Standardabweichung σO der Other-Scores SO
standardisiert werden, so verschiebt sich der Mittelwert der Other-Scores gegen 0, die
Standardabweichung der Other-Scores wird zu 1 und alle anderen Dichtefunktionen
skalieren entsprechend mit. Das Ergebnis sind auf Basis der Other-Scores standardi-
sierte Dichtefunktionen.
Die Abbildung 7.3 zeigt die Dichtefunktionen der Familie b.121.4.1 nach der Trans-
formation. Verglichen mit den vorherigen unstandardisierten Dichtefunktionen der Ab-
bildung 7.2 derselben Familie werden die Funktionen entsprechend skaliert. Der Mittel-
wert der Dichtefunktion fOther ist 0 und die Standardabweichung ist 1. Der Abstand der
Family-Dichtefunktion ist in Relation gleich geblieben, obwohl er nunmehr nur noch
ein Viertel des vorherigen Wertes hat. Gleiches gilt fur die Superfamily.
7. Bewertung der Ergebnisse 85
Abbildung 7.3: Standardisierte Dichtefunktionen des Reverse Corrected Scores fur dasProfil-MSA der Familie b.121.4.1
Mit dieser Standardisierung ist die Basis geschaffen, Scores uber Familien hinweg ver-
gleichen zu konnen. Vor allem aber besteht damit auch die Moglichkeit, die Veranderung
der Scoreverteilung zu beobachten, wenn Anderungen am Profil-MSA vorgenommen
werden. Die Frage die damit beantwortet werden soll lautet: Ist eine progressive Ex-
pansion des Profil-MSA sinnvoll?
7.4 Drift der Scores
Eine Anderung der HMM-Topologie durch Anderungen am Profil-MSA hat immer
auch eine Anderung der Bewertungsgrundlagen zur Folge, da sich innerhalb des Mo-
dells die Transitions- und Emissionswahrscheinlichkeiten verschieben. Durch die im Ab-
schnitt 6.2.4 eingefuhrte progressive Expansion des Profil-MSA kommt es zu derartigen
Anderungen des Profil-MSA in jedem Zyklus. Damit andern sich auch die Scoregroßen
und deren Verteilungen und es kommt zu einem Score-Drift3.
Beobachtet man die Dichtefunktionen nach jedem Zyklus, so stellt man je nach Aus-
gangsmodell, Anzahl und der Art der hinzugefugten Sequenzen fest, dass es teils leichte,
teils massive Verschiebungen in der Scoreverteilung gibt. Offensichtlich ist nicht jedes
Modell im gleichen Maße geeignet, Modellerweiterungen gewinnbringend aufzunehmen
3Der Begriff Drift wird der Nachrichten- oder Messtechnik entliehen. Dabei wird eineverhaltnismaßig langsame Anderung eines Wertes oder einer Eigenschaft eines Systems oder einerEinrichtung als Drift bezeichnet.
7. Bewertung der Ergebnisse 86
und damit die Erkennungsrate zu verbessern. Die Interpretation der Veranderungen an
den Dichtefunktionen soll der Beantwortung folgender Fragen dienen:
• Wie andert sich das Profil-HMM in Bezug auf die Scores?
• Sind mit den erweiterten Profil-HMMs bessere oder schlechtere Erkennungsraten
zu erwarten?
• Kommt es zu einer Uberanpassung des Modells oder zu einer Diversifikation des
Modells?
Die Anderung des Bewertungssystems zeigt sich am besten im direkten Vergleich der
Dichtefunktionen von einem Zyklus zum anderen. Wie im vorangegangenen Abschnitt
begrundet wurde, ist der Abstand der Dichtefunktionen der Family-, Superfamily- und
Other-Sequenzen ein wesentliches Merkmal in Bezug auf die erwarteten Erkennungsra-
ten. Liegen die Dichtefunktionen in Relation dicht beisammen, ist eine klare Zuordnung
schwierig; liegen sie weit auseinander, so ist eine klare Zuordnung einer Sequenz zu ei-
ner Family oder Superfamily einfacher. Ziel jedes Modells ist deshalb die Verwendung
eines Scores, der eine klare Zuordnung moglich macht und eine scharfe Trennung der
einzelnen Familien und Superfamilien vom Rest der Sequenzen ermoglicht. Ziel jeder
Veranderung am Modell ist aber auch, dass das Modell”zum Guten“ verandert wird,
so dass sich die Mittelwerte der Dichtefunktionen voneinander entfernen.
7.5 Graphische Bewertung der Scores
Die Abbildung 7.4 zeigt die Dichtefunktionen des Simple Corrected Score (SCS) und des
Reverse Corrected Score (RCS) nachdem die Sequenzen mit dem Profil-HMM der Fa-
milie e.3.1.1 einmal aligniert wurden. Die Abbildung zeigt, dass die Family-Sequenzen
der Familie e.3.1.1 uber die Scores eindeutig klassifiziert werden konnen und deshalb
hoch bewertet sind. Der standardisierte RCS liegt im Beispiel im Bereich von etwa 5 bis
30. Diese eindeutige Trennung der Family-Sequenzen vom Rest kann auch bei den ande-
ren Scores beobachtet werden. Die Sequenzen der Superfamily heben sich zwar sowohl
von den Family- also auch von den Other-Sequenzen ab, Uberlappungen auf beiden
Seiten machen eine eindeutige Klassifizierung der Sequenzen im Uberlappungsbereich
7. Bewertung der Ergebnisse 87
aber schwierig. Der RCS eignet sich in diesem Fall besser, als der SCS, da er - zumin-
dest im gegenstandlichen Testset - ein klareres Scoring in Abhangigkeit der Family und
Superfamily liefert.
Abbildung 7.4: Standardisierte Dichtefunktionen des Simple Corrected Scores und Re-verse Corrected Scores fur das Profil-MSA der Familie e.3.1.1 nach dem ersten Zyklus.
In zwei weiteren Zyklen wurde das Profil-MSA um je drei Sequenzen erweitert, wodurch
sich auch die Wahrscheinlichkeiten im HMM und damit die Bewertung der einzelnen
Sequenzen andern. Mit dem Profil-MSA wurden keine Sequenzen der gleichen Familie,
sondern Sequenzen anderer Familien aligniert, die trotzdem signifikante Ahnlichkeiten
mit dem Profil-HMM aufweisen. Der Grund dafur liegt in der verwendeten Testdaten-
bank, die keine weiteren Sequenzen der Familie e.3.1.1 mehr enthalten hat. Hinzu-
gefugt wurden stattdessen zwei Sequenzen der Familie e.3.1.2, zwei Sequenzen der
Familie e.3.1.3 und im dritten Zyklus je eine Sequenz der Familie d.165.1.1 und
c.100.1.1.
7. Bewertung der Ergebnisse 88
Abbildung 7.5: Standardisierte Dichtefunktionen des Simple und Reverse CorrectedScores fur das Profil-MSA der Familie e.3.1.1 nach dem dritten Zyklus.
Die Dichtefunktionen des SCS und RCS fur die Familie e.3.1.1 nach dem dritten
Zyklus zeigt die Abbildung 7.5. Deutlich wahrnehmbar ist die Verschiebung der Dich-
tefunktionen fur beide Scores aufgrund der Modellanderung. Die Dichtefunktionen der
Superfamily Sequenzen verschieben sich deutlich in Richtung Family, wodurch eine
deutlichere Trennung zu den Other-Sequenzen stattfindet, dem entgegen aber eine Un-
terscheidung zwischen Family und Superfamily nicht mehr moglich ist. Der Grund dafur
durfte sein, dass vier der hinzugefugten Sequenzen aus der gleichen Superfamily e.3.1
stammen wie das ursprungliche Profil-MSA und sich dadurch die Modellbeschreibung
mehr an der Superfamily orientiert. Nach drei Zyklen hat sich das Modell gegenuber
dem ursprunglichen Profil-HMM hinsichtlich der Erkennung von Superfamily Scores
deutlich verbessert. Sequenzen der gleichen Superfamily werden klar hoher bewertet
als Sequenzen abseits davon. Der Preis dafur ist eine geringere Trennscharfe zwischen
Family und Superfamily.
7. Bewertung der Ergebnisse 89
Weniger erfolgreich ist das Modell der Proteinfamilie a.138.1.1. Nach einem Zyklus
zeigt sich ein ahnliches Bild (Abbildung 7.6) wie im Beispiel zuvor. Die Sequenzen der
Family werden signifikant hoher bewertet als der Rest der Sequenzen. Die Sequenzen
mit der gleichen Superfamily heben sich von den Other-Sequenzen deutlich ab.
Abbildung 7.6: Standardisierte Dichtefunktionen des Reverse Corrected Scores fur dasProfil-MSA der Familie a.138.1.1 nach dem ersten Zyklus.
Wieder wurden in zwei Zyklen sechs weitere Sequenzen dem Profil-MSA angehangt
und wieder waren es drei Sequenzen der gleichen Superfamily und drei Sequenzen der
Familien g.49.1.1, g.44.1.1 und g.77.1.1. Das Modell entwickelte sich dadurch
aber nicht zu seinem Vorteil. Vergleicht man die Dichtefunktionen nach dem ersten
Zyklus mit denen nach dem dritten Zyklus (siehe Abbildung 7.7), so fallt auf, dass sich
die Dichtefunktion der Superfamily in Richtung Other bewegt und sich in der Form
ebenfalls der Verteilung der restlichen Sequenzen annahert.
Abbildung 7.7: Standardisierte Dichtefunktionen des Reverse Corrected Scores fur dasProfil-MSA der Familie a.138.1.1 nach dem dritten Zyklus.
7. Bewertung der Ergebnisse 90
Die Scores der Sequenzen der gleichen Familie werden kleiner, lassen aber im Vergleich
zu den Other-Sequenzen noch eine Klassifizierung zu. Die Standardabweichung der
Superfamily hat sich deutlich vergroßert, sodass sich die standardisierten Scores der
Superfamily nun uber den Bereich von -5 bis 8 erstrecken und eine Abgrenzung so-
wohl zu den Other- als auch zu den Family Sequenzen nicht uber alle Bereiche hinweg
moglich ist. Der Grund fur dieses Verhalten kann teilweise damit begrundet werden,
dass das ursprungliche Profil-MSA weniger Sequenzen hat und das Profil-HMM weniger
Emissionsspalten aufweist. In Relation dazu ist die Aufnahme von sechs entfernt ver-
wandten zusatzlichen Sequenzen in seinen Auswirkungen gewichtiger, als bei Modellen
mit vielen und langeren Sequenzen, zu denen einige kurze Sequenzen aligniert werden.
Ebenso ist eine teilweise”Unscharfe“ in der SCOP Klassifikation der Sequenzen der
Testdatenbank moglich. Es ist nicht klar, ob mit Hilfe von sequenzbasierten Methoden
die auf Struktur- und Funktionsinformation beruhende Klassifikation von SCOP (siehe
Kapitel 2.5.1) in allen Familien gleich gut abbildbar ist.
7.6 Numerische Bewertung der Scores
Die Beurteilung der multiplen Sequenzalignments uber alle Familien und alle Modelle
auf Basis von Dichtediagrammen ist aufwendig und nicht objektivierbar. Nach einigen
Testlaufen stellt sich heraus, dass bestimmte Modelle gut, andere wiederum schlechter
auf die progressive Expansion der Modelle reagieren. So spricht das Modell der Familie
b.6.1.1 und e.3.1.1 positiv auf das Hinzufugen neuer Sequenzen in den ersten zwei
Zyklen an, jenes von a.138.1.1 und a.3.1.4 jedoch nicht. Wie kann dieser Effekt in
Zahlen dargestellt werden?
Um diese Effekte numerisch auszudrucken, wird nach jedem Zyklus fur jede Family der
Mittelwert der standardisierten Family-, Superfamily- und Fold-Dichtefunktionen be-
rechnet und der Abstand in Relation zur Standardabweichung der Other-Dichtefunktion
gesetzt. Am Ende jedes Zyklus kann damit der Mittelwert uber alle Familien ermittelt
werden, um einen Gesamttrend abzulesen. Diese globalen Trendzahlen streichen zwar
keine einzelnen Modelle und Familien als besonders geeignet und ungeeignet heraus,
machen aber eine globale Tendenz uber alle Familien und Zyklen hinweg sichtbar.
7. Bewertung der Ergebnisse 91
Tabelle 7.2: Entwicklung der globalen Mittelwerte der Reverse Cor-rected Z-Scores und Simple Corrected Z-Scores der Families, Super-families und Folds in einem Testlauf uber vier Zyklen.
In der Tabelle 7.2 sind die Mittelwerte der Z-Scores aller Familien nach jedem von vier
Zyklen aufgelistet. Aus der Tabelle geht hervor, dass sich die Reverse Correded Z-Scores
(RCSZ) tendenziell verbessern. Die Families rucken bei jedem Zyklus zwar etwa naher
in Richtung Null, sind mit dem Abstand 18.09·σOther jedoch noch klar abgesetzt von den
Other-Sequenzen. Die Superfamilies rucken hingegen von der Other-Standardverteilung
(µ = 0, σ = 1) ab und vergroßern die Distanz im Laufe der Berechnungen auf 2.729.
Dass dies aber auch vom Scoretyp abhangig ist, zeigt der Vergleich mit dem Simple
Corrected Z-Scores (SCSZ). Im Laufe der vier Zyklen verandern sich die Family-Scores
ahnlich der des RCSZ , die Scores der Superfamily entwickeln sich aber schlechter. Die
Verteilungskurve der Superfamily driftet in die”falsche“ Richtung ab und der Abstand
der Mittelwerte verringert sich von 1.403 auf 1.318.
Die im Anhang C.3 gelisteten Auswertungen zeigen die Z-Scores auf Basis der Familie
am Ende des ersten und dritten Zyklus. Damit wird auch deutlich, wie unterschied-
lich und empfindlich die Modelle der Familien auf das Hinzufugen neuer Sequenzen
reagieren.
7.7 Streuung der Scores
In der im vorangegangen Abschnitt 7.4 dargestellten Untersuchung auf Basis der Score-
Dichtefunktionen wird ausschließlich auf einen einzigen Scoretyp Bezug genommen.
Dabei werden Scoretypen verwendet, deren Langenabhangigkeit gering ist und deren
Trennfahigkeit oder -scharfe sich als gunstig erwiesen haben. Nicht ohne Grund kommt
vermehrt der Reverse Corrected Score zur Anwendung, wie Lackner et al. schon in
7. Bewertung der Ergebnisse 92
[27] belegen. Wie die Auswertungen der Scores und deren Entwicklung uber mehrere
Zyklen aber auch zeigen, ist die Trennscharfe der verwendeten Scores sowohl durch
den Scoretyp als auch durch das Profil-MSA begrenzt. Die Sequenzen gleicher Fami-
lies setzen sich meist deutlich vom Rest der Sequenzen ab, aber schon die eindeutige
Klassifikation der Superfamily ist nicht immer gewahrleistet.
Eine Gegenuberstellung der Corrected Scores zeigen Streudiagramme in der Abbildung
7.8, in denen auf den X- und Y-Achsen jeweils die Werte eines Scoretyps aufgetragen
werden. Dargestellt werden unterschiedliche Scorepaarungen eines Modells der Familie
a.1.1.2. Die cyanfarbigen und roten Punkte zeigen die Scores der Superfamilies und
Families und die blauen Punkte die Other-Scores. Mit einem farbigen Kreuz sind jene
Scores markiert, die das Profil-MSA beim Start der Zyklen definiert haben, und die
grunen Punkte deuten jene Scores an, die im Laufe der sechs Zyklen dem Profil-MSA
hinzugefugt wurden.
Abbildung 7.8: Ausgewahlte Scores der Familie a.1.1.2 in Streudiagrammen.
Im Vergleich der Diagramme zeigt sich, dass einige Scores nahezu identische Bewer-
tungen liefern, wie beispielsweise die Kombinationen aus Simple Corrected Score und
HMM Corrected Score. Die Scores fadeln sich mit geringer Streuung entlang eines engen
Kanals auf. Andere Diagramme zeigen zwar deutlichere Streuungen, trotzdem aber im
Verlauf noch deutliche Abhangigkeiten. In allen Scorekombinationen heben sich die Fa-
milies gut und die Superfamilies teilweise vom Rest der Sequenzen ab, was aufgrund der
Erfahrungen aus den Dichtediagrammen keine Uberraschung darstellt. Eine komplette
7. Bewertung der Ergebnisse 93
Zusammenstellung aller Scorepaarungen der Zyklen 1 und 3 fur die Familie a.1.1.2
zeigen die Abbildungen A.2 und A.3 im Anhang A4.
Aus den Diagrammen wird deutlich, dass keine Kombination aus zwei Scores die
Trennscharfe wesentlich verbessern wurde. Dies bestatigt auch das folgende Diagramm.
Abbildung 7.9: Transformierte Scores der Familie a.1.1.2 nach einer Hauptkomponen-tenanalyse nach dem dritten aus insgesamt sechs Zyklen.
Die Abbildung 7.9 zeigt ein Streudiagramm mit dem Ergebnis einer Hauptkomponen-
tenanalyse (Principle Components Analysis ; PCA), ein Verfahren der multivariaten
Statistik. Ziel der PCA ist - stark vereinfacht und auf diesen Anwendungsfall bezogen
erklart - mehrdimensionale, teils stark korrelierende Datensatze so zu komprimieren,
dass die Datensatze in moglichst wenigen Dimensionen und unkorreliert dargestellt
werden konnen und moglichst viel der ursprunglichen Information erhalten bleibt. Im
konkreten Fall wurden die 6-dimensionalen Scoredaten auf 2-Dimensionen reduziert,
um sie im Streudiagramm visualisieren zu konnen.
4Eine Zusammenstellung von Streudiagrammen aller Familien uber 6 Zyklen befindet sich auf dembeigelegten Datentrager (siehe Anhang E)
7. Bewertung der Ergebnisse 94
Das Diagramm in der Abbildung 7.9 belegt die Abhangigkeiten der Scores auch, ver-
gleicht man es beispielsweise mit den Reverse Corrected Score Diagrammen aus Abbil-
dung 7.8. Wie aus dem Diagramm deutlich wird, unterscheidet sich das PCA Diagramm
kaum vom Ergebnis anderer Scorepaarungen. Paarungen mit dem Reverse Corrected
Score und dem Reverse Correted Score (Foreward) zeigen durchaus vergleichbare Er-
gebnisse. Dies ist insofern bemerkenswert, als die PCA eine Linearkombination aller
sechs Dimensionen bildet. Dies bestatigt aber auch die Annahme, dass die Abhangigkeit
der Scores untereinander groß ist und schon ein bis zwei Scores die Informationen na-
hezu vollstandig enthalten5.
Die Voraussetzungen fur mindestens zwei elementare unabhangige Merkmale ist damit
kaum gegeben. Keine der dargestellten Scorepaarungen wurde demnach in der unter-
suchten Familie wesentlich zur Bildung eines trennfahigen nichtelementaren Merkmals
geeignet sein.
7.8 Entropie der Emissionsmatrix
Zur Uberprufung von Scoring-Matrizen wird vielfach der Informationsgehalt oder die
Informationsdichte der Matrizen festgestellt. Diese gibt - vereinfacht ausgedruckt -
daruber Auskunft, inwieweit eine Scoring-Matrix Informationen uber die Verteilung
von Aminosauren in Sequenzen widerspiegelt oder ob eine Matrix eine reine Zufalls-
verteilung beschreibt, also keinerlei verwertbare Information innehat. Grundlage der
Bemessung des Informationsgehalts ist die Annahme, dass die Information, die ei-
ne Nachricht enthalt, umgekehrt proportional zu ihrer Eintrittswahrscheinlichkeit ist.
Tritt innerhalb einer Sequenz ein Zeichen mit der Wahrscheinlichkeit p = 1.0 auf, so
enthalt dieses Ereignis keinerlei Information. Ist ein Ereignis unwahrscheinlich und tritt
es doch ein, so vermittelt dies eine Information großer Null. Die gleiche Annahme gilt
bei Scoring-Matrizen fur die einzelnen Zeichenpaarungen, die darin abgebildet werden.
5Genauere Vergleiche der PCA Ergebnisse zeigen, dass mit den meisten getesteten HMMs in einemScore schon 90% und mit zwei Scores 98% der Varianz beschrieben werden. Werden im Laufe derIterationen Sequenzen dem Profil-MSA hinzugefugt, so sinkt dieser Wert nach 6 Zyklen auf etwa 80%mit einem Score und 95% mit zwei Scores.
7. Bewertung der Ergebnisse 95
Gemessen wird die Informationsdichte mit der von C.E. Shannon6 eingefuhrten Entro-
pie, einem quantitativen Maß zur Bemessung des mittleren Informationsgehalts eines
Zeichensystems. Ist die Entropie groß, so ist der Informationsgehalt des Systems groß,
geht die Entropie gegen Null, so geht auch der Informationsgehalt gegen Null.
Wird fur die Berechnung der mittleren Entropie H der Logarithmus mit der Basis 2
verwendet, dann ist die Einheit der mittleren Entropie H mit bit/Symbol definiert. Die
mittlere Entropie H als das Maß einer Wahrscheinlichkeitsverteilung P = (p1, . . . , pn)
selbst ist definiert als [29]:
H(P ) := −n∑
i=1
pi · log2(pi) in bit/Symbol (7.2)
Der Informationsgehalt I eines einzelnen Symbols mit dem Logarithmus der Basis 2
gerechnet ergibt die Dimension bit und ist definiert als:
Ii := −log2(pi) in bit (7.3)
Da alle pi immer kleiner 1 sind, ergibt sich aus obiger Formel, dass die negative Summe
der negativen Logarithmuswerte immer einen positiven Wert fur H bilden.
Hermando et al. beschreiben in [20], wie die Entropie eines Hidden Markov Modells
exakt und mit erheblichen Aufwand berechnet werden kann. Eine in dieser Arbeit stark
vereinfachte Methode berucksichtigt ausschließlich die Entropie der Emissionsmatrix,
um die Anderungen der Entropie im Laufe mehrerer Zyklen zu dokumentieren. Ziel
ist die Feststellung, ob die Entropie durch die Aufnahme neuer Sequenzen ab- oder
zunimmt und ob ein Zusammenhang der Anderungen der Wahrscheinlichkeitsdichten
mit den Entropieanderungen hergestellt werden kann.
Nimmt man die Emissionsmatrix eines Profil-HMM, so gibt jede Emissionsspalte der
6Claude Elwood Shannon, 1916 - 2001, US-amerikanischer Elektroingenieur und Mathematiker; Be-grunder der Informationstheorie und wesentliche Beitrage zur Kryptographie und Nachrichtentechnik;er verfasste 1948 die Mathematical Theory of Communication und 1949 die Communication Theoryof Secrecy Systems [1].
7. Bewertung der Ergebnisse 96
Matrix die Wahrscheinlichkeitsverteilungen der 20 Zeichen des Alphabets fur Ami-
nosauren ΣA wieder. Die Summe der Eintrittswahrscheinlichkeiten aller einzelnen Zei-
chen einer Emission ergibt demnach immer 1. HMModeler lasst auch Summen großer
1 zu [39], diese Besonderheiten kommen hier aber nicht zur Anwendung. Folgt man
obiger Definition der Entropie, so lasst sich fur jede Spalte der Informationsgehalt be-
rechnen und fur die gesamte Matrix der Mittelwert der Spaltenentropien ermitteln.
Diese Große kann als die Entropie der Emissionsmatrix angesehen werden.
Wird im Laufe der Berechnungen in jedem Zyklus das Profil-HMM um einige Sequenzen
erweitert, so andern sich dadurch die Emissions- und Transitionswahrscheinlichkeiten
des Modells ebenso, wie die Alignierungseigenschaften, die Dichtefunktionen und die
Entropie der Emissionsmatrix. Die Tabelle 7.3 am Ende dieses Abschnitts stellt die
Anderung der Z-Scores und die Anderung der Entropie zwischen den Zyklen 1 und 3
gegenuber. Die grun markierten Z-Scores der Superfamily weisen auf eine Vergroßerung
(=Verbesserung) der Mittelwertabstande hin, rot markierte Z-Scores der Superfamily
weisen auf eine Verringerung (=Verschlechterung) der Mittelwertabstande hin. Stellt
man diese Bewegungen den Trends der Entropie gegenuber, so zeigt sich, dass die
Entropie in vielen Fallen zwar steigt, sich die Steigerung des Informationsgehalts aber
nicht immer positiv auf den Superfamily Z-Score auswirkt. Der Grund dafur durfte
sein, dass die Aufnahme neuer Sequenzen das HMM zwar mit Information anreichert,
diese Informationen jedoch nicht ausschließlich auf die Family oder Superfamily bezo-
gen sind. Werden dem Modell neue, im Verwandtschaftsverhaltnis entfernter liegende
Sequenzen hinzugefugt, so diversifiziert das Modell, wie auch in der Abbildung 7.10
deutlich wird, indem die Streuung der Scores steigt. Werden Sequenzen der gleichen
Familie hinzugefugt, so spezialisiert sich das Modell auf diese eine Familie.
Auffallend ist die Steigerung der Entropie der Familien a.1.1.2 und a.138.1.1 (siehe
Tabelle 7.3, Spalte deltaE). Im Fall der Familie a.1.1.2 wurden die Thresholds zur
Sequenzvorauswahl wirksam, so dass in den drei Zyklen nur zwei neue Sequenzen dem
Modell hinzugefugt wurden, wovon eine Sequenz aus der gleichen Superfamilie stammt.
Im Fall der Familie a.138.1.1 wurden vier neue Sequenzen hinzugefugt, drei davon
stammen aber aus der gleichen Familie. Die Entropie steigt damit deutlich, das heißt,
das HMM hat seitens der Emissionsmatrix an Information gewonnen. Trotzdem fuhrt
7. Bewertung der Ergebnisse 97
(a) Zyklus 1 (b) Zyklus 6
Abbildung 7.10: Anderung der PCA-transformierten Scores der Familie d.32.1.2 infolgeder Expansion des Profil-MSA zwischen dem ersten und dem sechsten Zyklus.
diese Steigerung zu keinem merkbaren Trend in den Z-Scores. Im ersten Fall verbessert
sich das Ergebnis minimal, im zweiten verschlechtert sich das Ergebnis sogar. Der
Grund dafur durfte sein, dass im Zyklus 3 in der Testdatenbank keine neuen Sequenzen
zur Verfugung standen, auf die das Modell abgestimmt war. Alle zum Modell passenden
Sequenzen waren zu diesem Zeitpunkt schon Teil des Modells (siehe Tabelle 7.1).
7. Bewertung der Ergebnisse 98
Tab
elle
7.3:
Ver
glei
chder
rela
tive
nA
nder
unge
nder
Dic
hte
funkti
onen
mit
den
Entr
opie
nzw
isch
enden
Zykle
n1
und
3.
7. Bewertung der Ergebnisse 99
7.9 Zusammenfassung der Evaluationsergebnisse
Die Parametrisierung der HMMs als Basis einer objektiven und vergleichbaren Be-
wertung der Alignments kann als kritisch und schwierig erachtet werden. Die bis dato
ublichen Scores, die nur Singlealignments mit dem Profil-MSA bewerten, sind dafur nur
”begrenzt verwertbar“. Einige davon sind sowohl modell- als auch sequenzlangenab-
hangig, womit eine modellubergreifende Bewertung und ein objektiver Vergleich der
Alignierungsergebnisse kaum moglich ist. Die besten Ergebnisse konnen mit dem Re-
verse Corrected Scores erreicht werden, die sich als die stabilsten Scores herausstellten.
Keine Kombination aus zwei Scores stellt jedoch eine wesentliche Verbesserung der
Trennscharfe in Aussicht.
Uber die Verwendung von langen(teil-)korrigierten Z-Scores und die Standardisierung
der Dichtefunktionen, konnen einige Nachteile der Scores abgeschwacht werden. Ver-
gleicht man die Ergebnisse der Alignierungen uber mehrere Zyklen hinweg so fallt aber
auf, dass der Erfolg oder Misserfolg auch wesentlich vom Profil-MSA abhangt. Ein
HMM, das mit einem Profil-MSA mit wenigen Sequenzen beschrieben wird, reagiert
sensibler auf das Hinzufugen familienfremder Sequenzen als ein Modell, das mit einem
umfangreichen MSA trainiert wird. Damit steht auch fest, dass die Expansion eines
Modells nur dann sinnvoll ist, wenn das Modell stabil genug ist, um Sequenzen ent-
fernter Verwandter aufzunehmen. Ist dies nicht der Fall, so driftet das Modell mitunter
so weit ab, dass es eine Familie oder einen”Clan“7 nur noch unzureichend beschreibt.
Ein direkter Zusammenhang der Entropie der Emissionsmatrizen mit den Scoringer-
gebnissen kann nicht festgestellt werden. Die Anderung der Entropie kann mitunter
als Zeichen der Diversifikation oder Spezialisierung eines HMM interpretiert werden,
direkte Auswirkungen auf die Alignierungsergebnisse konnen mit den aktuellen Test-
daten aber nicht nachgewiesen werden. Ein Grund fur teilweise stark unterschiedliche
Testergebnisse ist mitunter in den Testdaten selbst zu suchen, die bestimmte Modelle
aufgrund der Zusammensetzung der Daten begunstigen oder andere benachteiligen.
7Der Begriff des Clans wird hierbei in Bezug auf die Datenbank Pfam verwendet.”A recent deve-
lopment in Pfam has enabled the grouping of related families into clans.“[5]
8
Zusammenfassung und Ausblick
Eines der wesentlichen Anliegen der Bioinformatik ist die effiziente Berechnung der
Ahnlichkeit von Sequenzen. Die Losung dieser Aufgabe ist insofern von hoher Relevanz,
als aus der Ahnlichkeit von Aminosauresequenzen auf die evolutionare Verwandtschaft
von Proteinen und daruber hinaus auf deren strukturelle und funktionelle Ahnlichkeit
geschlossen werden kann. Multiple Sequenzalignments eignen sich besonders, die in den
Sequenzen vorkommenden Muster und Ahnlichkeiten zu identifizieren. Deren Erstel-
lung gilt aber als nicht trivial und algorithmisch komplex.
Im Rahmen dieser Arbeit wurde eine Methode entwickelt, mit der auf Basis von Hidden
Markov Modellen multiple Sequenzalignments erzeugt werden. Die folgende Zusam-
menfassung beschreibt die Implementierung, mit der diese Aufgabe gelost wird und
die Ruckschlusse, die aus den Testergebnissen gezogen werden. Den Abschluss bildet
der Ausblick mit Vorschlagen, wie die Ergebnisse in Zukunft weiter verbessert werden
konnten.
8.1 Zusammenfassung
Die vorliegende Thesis zeigt die Anwendung und Evaluation eines Profil Hidden Mar-
kov Modells (Profil-HMM) zur Berechnung eines multiplen Sequenzalignments. Dabei
wird eine mehrstufige Strategie implementiert, mit der die Sequenzen nach der Einzela-
lignierung mit dem Profil-MSA in einem iterativen Prozess zu einem Vereinigungs-MSA
100
8. Zusammenfassung und Ausblick 101
(VMSA) aligniert werden. In einem zusatzlichen Arbeitsgang wird mit einer alternati-
ven Alignierungsmethode das VMSA nachaligniert, um das Ergebnis weiter zu verbes-
sern. Dazu werden die vom HMM unaligniert gebliebenen Bereiche in Snippets extra-
hiert und diese mit ClustalW aligniert. Danach werden die alignierten MSA-Snippets
wieder in das VMSA zuruckgefuhrt.
Eine Erweiterung erfahrt die Methode, indem das Profil-MSA in einem iterativen Pro-
zess schrittweise mit Sequenzen aligniert und erweitert wird, um es in jedem Durchlauf
neuerlich fur Alignierungen wiederzuverwenden.
Die Trennfahigkeit, Sensibilitat und Stabilitat eines Hidden Markov Modells und die
Alignierungsergebnisse hangen von einer Reihe von Parametern ab. Die Sequenzanzahl
und die Sequenzlange der Profil-MSA spielt dabei eine ebenso wesentliche Rolle, wie
die Anzahl der Emissionsspalten und die Sequenzlange der zu alignierenden Sequen-
zen. Diese Multivariabilitat der Parameter erschwert einerseits die Vorhersagbarkeit
der Modelleigenschaften, anderseits auch die Vergleichbarkeit der Testergebnisse uber
mehrere Modelle und Zyklen hinweg.
Da die quantitative Beurteilung eines MSA aufgrund fehlender Gesamtscores schwer
moglich ist, stutzt sich die Bewertung der Ergebnisse auf die Beobachtung der Sco-
res und deren dynamische Entwicklung im iterativen Prozess. Dabei wird - vorerst
graphisch - die dynamische Entwicklung der Scoredichtefunktionen wahrend mehre-
rer Iterationen bewertet. Eine nachfolgende numerische Bewertung erfolgt uber den
Vergleich der Entropieanderung der Emissionsmatrizen und der Mittelwertanderung
standardisierter Z-Scores.
Die Implementation und Evaluation der Ergebnisse zeigt, dass die Erstellung eines
multiplen Sequenzalignments uber ein Hidden Markov Modell grundsatzlich und pro-
blemlos moglich ist. Kritisch ist vielmehr die Parametrisierung eines HMM in Form
des Trainings und die objektive Bewertung der Ergebnisse. Die zur Verfugung stehen-
den Scores, die derzeit nur Singlealignments mit dem HMM bewerten, sind dafur nur
begrenzt verwertbar. Diese sind teilweise stark modell- und sequenzlangenabhangig,
womit eine modellubergreifende Bewertung und ein objektiver Vergleich der Alignie-
rungsergebnisse schwer moglich ist.
8. Zusammenfassung und Ausblick 102
8.2 Ausblick
Diese Arbeit belegt einerseits die Moglichkeiten der Anwendung eines HMMs im Se-
quenzalignment, andererseits aber auch die Grenzen der derzeitigen Implementierung.
Es empfiehlt sich, in der Fortsetzung dieser Arbeit eine Reihe von Verbesserungen im
Bereich der Algorithmik, dem Scoring, der Erstellung der Profil-MSA, sowie in die
Performance einfließen zu lassen.
Wie die Streudiagramme und die PCA mit den Beispieldaten belegen, stellt keine
Kombinationen der Scores der gepruften Familien eine wesentliche Verbesserung der
Trennscharfe in Aussicht. Die Analyse der PCA belegt aber auch, dass je spezialisierter
ein Modell ist, desto eher ein einziger Scores die Gesamtvarianz des Systems erklart,
und je diversifizierter ein Modell ist, desto eher sich die Varianz auf mehrere Scores
verteilt. Demnach ist durchaus denkbar, dass bei bestimmten Modellen und bestimmten
Familien die Kombination von zwei Scores Vorteile bringen konnte.
Eine Untersuchung der Testdaten und Profil-MSA zeigt, dass die Profil-MSA aus-
schließlich aus allen Sequenzen einer einzigen Familie beschrieben werden, die sich auch
in den Testdaten wiederfinden. Das Problem liegt darin, dass dieselben Sequenzen bei
der Klassifikation und beim Test des Verfahrens zur Anwendung kommen, die schon
bei der Profil-MSA-Bildung beteiligt waren. Um dieses Defizit zu beheben, konnte die
sogenannte Leave-one-out Methode angewendet werden, bei dem jeweils jene Sequenz
aus dem Profil-MSA isoliert wird, die mit dem verbleibenden MSA aligniert wird.
Nachdem die Erfolgsaussichten einer guten Alignierung wesentlich von der Anpassung
der Parameter abhangt und dafur oft mehrere Testlaufe notwendig sind, ist fur einen
praktischen Einsatz eine Verbesserung der Performance wunschenswert. Derzeit sind
Rechenzeiten von einigen Stunden bis Tagen durchaus ublich. Eine Parallelisierung der
Software oder die Umstellung auf ein verteiltes Softwaresystem konnte die Performan-
ce deutlich erhohen. Eine Verteilung der Prozesse in einem Cluster ware mit einem
vergleichsweise geringen Aufwand implementierbar, ohne große Anderungen an den
Applikationen selbst vornehmen zu mussen.
Literaturverzeichnis
[1] Encyclopedia Britannica - Online Ausgabe, 2010. URL: http://www.britannica.
com/EBchecked/topic/365793/Andrey-Andreyevich-Markov (Stand: 28. Juli
2010)).
[2] A.G. Murzin et al.: SCOP: A Structural Classification of Proteins Database for the
Investigation of Sequences and Structures. Journal of Molecular Biology, 247:536-
540, 1995.
[3] G.R. Cochrane et al.: Petabyte-scale innovations at the European Nucleotide
Archive. In Nucleic Acids Research. Oxford University Press, Oxford, October
2008.
[4] L. Lo Conte et al.: SCOP database in 2002: refinements accommodate structural
genomics. Nucleic Acids Research, 230:264–267, 2002.
[5] R.D. Finn et al.: Pfam: clans, web tools and services. Nucleic Acids Res, 34:247–
251, 2006.
[6] F. Auer: Scoring Schemes and Parameter Prediction for Profile HMMs. Diplo-
marbeit, Fachhochschule Salzburg; Fachbereich Informationstechnik und System-
Management, Puch/Salzburg, 2009.
[7] D. Bindreither, S. Wegenkittl, F. Auer, and P. Lackner: Expert knowledge enhanced
structure based profile HMMs for protein sequence families. In Posterpresentation
at the German Conference on Bioinformatics 2009, 2009.
[8] L. Bonetta: Genome sequencing in the fast lane. In Nature Methods, volume 3,
pages 141–147. New York, February 2006.
103
Literaturverzeichnis 104
[9] G.R. Cochrane and M.Y. Galperin: The 2010 Nucleic Acids Research Database
Issue and online Database Collection: a community of data resources, volume 38.
Oxford University Press, January 2010.
[10] M.O. Dayhoff and R.M. Schwartz: Matrices for detecting distant relationships. In
Atlas of protein sequence and structure, volume 5, pages 353–358, Washington
D.C., 1978. National Biomedical Research Foundation.
[11] M.O. Dayhoff, R.M. Schwartz, and B.C. Orcutt: A model of evolutionary change
in proteins. In Atlas of Protein Sequence and Structure, volume 5, pages 345–352,
Washington D.C., 1978. National Biomedical Research Foundation.
[12] R. Durbin, S. Eddy, A. Krogh, and G. Mitchison: Biological sequence analysis:
probabilistic models of proteins and nucleic acids. Cambridge University Press,
Cambridge, 11th edition, 2006.
[13] S. Eddy: Hidden Markov Models And Large-Scale Genome Analysis. 1997.
[14] D.F. Feng and R.F. Doolittle: Progressive Sequence Alignment as a Prerequisite
to Correct Phylogenetic Trees. Journal of Molecular Evolution, 25:351-360, 1987.
[15] G.A. Fink: Mustererkennung mit Markov-Modellen. Teubner B.G. GmbH, Stutt-
gart, 2003.
[16] W. Gilbert and A.M. Maxam: A new method for sequencing DNA. In Proceedings
of the National Academy of Sciences U.S.A., volume 74, pages 560–564, Washing-
ton, DC, February 1977.
[17] D. Gusfield: Algorithms on Strings, Trees, and Sequences. Cambridge University
Press, Cambridge, 1997.
[18] D. Haussler, A. Krogh, I.S. Mian, and K. Sjolander: Protein Modeling using Hidden
Markov Models: Analysis of Globins. Technical report, University of California at
Santa Cruz, Santa Cruz, USA, 1993.
[19] S. Henikoff and J.G. Henikoff: Amino acid substitution matrices from protein
blocks. In Proceedings of the National Academy of Sciences of the United States
Literaturverzeichnis 105
of America, volume 89, pages 10915–10919, Washington, DC, 1992. Proceedings
of the National Academy of Sciences, USA.
[20] D. Hernando, V. Crespi, and G. Cybenko: Efficient Computation of the Hidden
Markov Model Entropy for a Given Observation Sequence. In IEEE Transaction
on Information Theory, volume 51, pages 2681–2685. IEEE, July 2005.
[21] M. T. Hutt und M. Dehnert: Methoden der Bioinformatik: Eine Einfuhrung. Sprin-
ger Verlag, Berlin, 2006.
[22] IUPAC-IUBMB: Joint Commission on Biochemical Nomenclature and Nomen-
clature Commission of IUBMB: Nomenclature and Symbolism for Amino Acids
ans Peptides Biochemical Nomenclature and Related Documents. International
Union of Pure and Applied Chemistry and International Union of Biochemistry
and Molecular Biology / Joint Commission on Biochemical Nomenclature, UK,
1984. URL: http://www.iupac.org/publications/pac/1984/pdf/5605x0595.
pdf (Stand: 23.Marz 2010).
[23] M.Y. Kao: Encyclopedia of Algorithms. Springer Verlag, New York, 1st edition,
2008.
[24] S. Karlin and S.F. Altschul: Methods for assessing the statistical significance of
molecular sequence features by using general scoring schemes. Proceedings of the
National Academy of Sciences, USA, 87:2264-2268, 1990.
[25] G. Kramer: Kleines Lexikon der Epileptologie (Taschenbuch). Thieme Verlag,
Stuttgart, 1. Auflage, 2005.
[26] A. Krogh, M. Brown, M. Mian, M. Sjolander, and D. Haussler: Hidden Markov
models of computational biology: Applications to protein modelling. Journal of
Molecular Biology, 235:1501-1531, 1994.
[27] P. Lackner, F. Auer, M. Radlingmaier, and S. Wegenkittl: Optimierte Modelle
zur Beschreibung von Proteinfamilien. In Proceedings der FFH2009, Drittes
Forschungsforum der Osterreichischen Fachhochschulen, April 15–16, 2009. Fach-
hochschule Karnten, 2009.
Literaturverzeichnis 106
[28] D.J. Lipman and W.R. Pearson: Rapid and sensitive protein similarity searches.
Science, 227:1435-1441, 1985.
[29] D. Lochmann: Digitale Nachrichtentechnik: Signale, Codierung,
Ubertragungssysteme, Netze). Verlag Technik, Berlin, 2. Auflage, 1997.
[30] C.D. Manning, P. Raghavan, and H. Schutze: An Introduction to Informa-
tion Retrieval. Cambridge University Press, Cambridge, 2009. Online Aus-
gabe URL: http://nlp.stanford.edu/IR-book/pdf/irbookonlinereading.
pdf (Stand: 03.Marz 2010.
[31] K.J. Menlove, M. Clement, and K.A. Crandall: Similarity Searching Using BLAST.
Methods in Molecular Biology. Humana Press / Springer, Provo, 2009.
[32] R. Merkl und S. Waack: Bioinformatik Interaktiv: Algorithmen und Praxis. Wiley-
VCH Verlag GmbH, Weinheim, 2. Auflage, 2009.
[33] D.W. Mount: Bioinformatics: Sequence and Genome Analysis. Cold Spring Harbor
Laboratority Press, New York, 2001.
[34] NCBI - National Center for Biotechnology Information: FormatDB and FastaCmd,
2007. URL: http://www.ncbi.nlm.nih.gov/staff/tao/URLAPI/formatdb_
fastacmd.html#t1.1 (Stand: 25. Juli 2010).
[35] S.B. Needleman and C.D. Wunsch: A general method applicable to the search for
similarities in the amino acid sequence of two proteins. Journal of Molecular
Biology, 48:443-453, 1970.
[36] Nobelprize.org - The Official Web Site of the Nobel Prize: The Nobel Prize in
Chemistry 1980 - Paul Berg, Walter Gilbert, Frederick Sanger, 1980. URL:
http://nobelprize.org/nobel_prizes/chemistry/laureates/1980/ (Stand:
25. Juli 2010).
[37] L. Papula: Mathematik fur Ingenieure und Naturwissenschaftler, Band 3. Viewweg
Verlag, Braunschweig, 4. Auflage, 2001.
[38] W. Pirovano and J. Heringa: Multiple Sequence Alignment, volume 452 of Methods
in Molecular Biology. Humana Press / Springer, Totowa, NJ, May 2008.
Literaturverzeichnis 107
[39] M. Radlingmaier: Optimierte Hidden Markov Modelle fur das Sequenz-Scoring auf
Basis von Multiplen Alignments. Diplomarbeit, Fachhochschule Salzburg; Fach-
bereich Informationstechnik und System-Management, Puch/Salzburg, 2007.
[40] F. Sanger, S. Nicklen, and A.R. Coulsen: DNA sequencing with chain-terminating
inhibitors. In Proceedings of the National Academy of Sciences U.S.A., volume
74, pages 5463–5467, Washington, DC, October 1977.
[41] R. Sedgewick: Algorithmen. Pearson Education Inc., Munchen, 2. Auflage, 2002.
[42] K. Sharma: Bioinformatics: Sequence Alignment and Markov Models. McGraw-
Hill Professional, New York, 1st edition, 2008.
[43] T.F. Smith and M.S. Waterman: Identification of Common Molecular Subse-
quences. Journal of Molecular Biology, 147:195-197, 1981.
[44] E.L. Sonnhammer, S.R. Eddy, and R. Durbin: Pfam: A Comprehensive Database
of Protein Domain Families Based on Seed Alignments. PROTEINS: Structure,
Function, and Genetics, 28:405-420, 1997.
[45] J.D. Thompson, D.G Higgins, and T.J. Gibson: CLUSTAL W: improving the
sensitivity of progressive multiple sequence alignment through sequence weighting,
positionspecific gap penalties and weight matrix choice. Nucleic Acids Research,
22:4673-4680, 1994.
[46] C.S. Tsai: An Introduction to Computational Biochemistry. Wiley-Liss, Inc, New
York, 2002.
[47] S. Wegenkittl, F. Auer, D. Bindreither, and P. Lackner: Expert knowledge enhanced
structure based profile HMMs for protein sequence families. In Posterpresentation
at the 3DSig 2009, 2009.
Abkurzungsverzeichnis
BLOSUM . . . . . . . . . . . . . . Block Substitution Matrix
CCS . . . . . . . . . . . . . . . . . . . Concise Classification String
CDF . . . . . . . . . . . . . . . . . . . Cumulated Distribution Function
DNA . . . . . . . . . . . . . . . . . . Desoxyribonukleinsaure
FASTA . . . . . . . . . . . . . . . . Fasta Dateiformat
HCS . . . . . . . . . . . . . . . . . . . HMMer Corrected Score
HCSF . . . . . . . . . . . . . . . . . HMMer Corrected Score (Forward)
HMM . . . . . . . . . . . . . . . . . . Hidden Markov Modell
MSA . . . . . . . . . . . . . . . . . . Multiple Sequence Alignment
NCDF . . . . . . . . . . . . . . . . . NormCDF based on Mean
NCDFM . . . . . . . . . . . . . . . NormCDF based on Median
PAM . . . . . . . . . . . . . . . . . . Point Accepted Mutations
PCA . . . . . . . . . . . . . . . . . . . Principle Components Analysis
PDB . . . . . . . . . . . . . . . . . . . Protein Data Bank
PFAM . . . . . . . . . . . . . . . . . Protein Families
PHMM . . . . . . . . . . . . . . . . Profile Hidden Markov Modell
PIG . . . . . . . . . . . . . . . . . . . Protein Identification Group
PSC . . . . . . . . . . . . . . . . . . . Plain Score
PSA . . . . . . . . . . . . . . . . . . . Pairwise Sequence Alignment
PSF . . . . . . . . . . . . . . . . . . . Plain Score (Forward)
REGEX . . . . . . . . . . . . . . . Regular Expression
RCS . . . . . . . . . . . . . . . . . . . Reverse Corrected Score
RCSF . . . . . . . . . . . . . . . . . . Reverse Corrected Score
RNA . . . . . . . . . . . . . . . . . . Ribonukleinsaure
108
Abkurzungsverzeichnis 109
RSC . . . . . . . . . . . . . . . . . . . Reverse Score
RSF . . . . . . . . . . . . . . . . . . . Reverse Score (Forward)
SCCS . . . . . . . . . . . . . . . . . . Set of Concise Classificaton Strings
SCOP . . . . . . . . . . . . . . . . . Structural Classification Of Proteins
SCS . . . . . . . . . . . . . . . . . . . Simple Corrected Score
SCSF . . . . . . . . . . . . . . . . . . Simple Corrected Score (Forward)
SQS . . . . . . . . . . . . . . . . . . . Sequenzstring
SQN . . . . . . . . . . . . . . . . . . . Sequenzname
SQL . . . . . . . . . . . . . . . . . . . Structured Query Language
SSC . . . . . . . . . . . . . . . . . . . Simple Score
SSA . . . . . . . . . . . . . . . . . . . Single Sequence Alignment
SSF . . . . . . . . . . . . . . . . . . . Simple Score (Forward)
VMSA . . . . . . . . . . . . . . . . . Vereinigungs-MSA
Anhang
110
A
111
A. Tabellen und Abbildungen 112
Tabellen und Abbildungen
A.1 PAM250 Matrix
Tab
elle
A.1
:D
ieP
AM
250
Mat
rix
zeig
tdie
Wah
rsch
einlich
keit
enei
ner
durc
hei
ne
Muta
tion
ents
tanden
enD
ista
nz
von
250
PA
M-E
inhei
ten
inP
roze
nt
(bas
iere
nd
auf
[11]
).D
ieT
abel
leze
igt
bei
spie
lsw
eise
,das
sb
eim
Ver
glei
chzw
eier
Am
inos
aure
sequen
zen,
ander
ener
ster
Pos
itio
nzu
vor
ein
A(A
la)
stan
d,
mit
einer
Wah
rsch
einlich
keit
von
13%
nac
hei
ner
akze
pti
erte
nM
uta
tion
eben
sonoch
ein
A(A
la)
steh
enw
ird.
Eb
enso
bes
teht
eine
3%ig
eC
han
ce,
das
sb
eiei
ner
Muta
tion
das
Azu
einem
R(A
rg)
muti
ert.
A. Tabellen und Abbildungen 113
A.2 Score-Streudiagramme
Tab
elle
A.2
:Str
eudia
gram
me
der
Sco
res
der
Fam
ilie
a.1.
1.2
nac
hdem
erst
enZ
yklu
s.
A. Tabellen und Abbildungen 114
Tab
elle
A.3
:Str
eudia
gram
me
der
Sco
res
der
Fam
ilie
a.1.
1.2
nac
hdem
dri
tten
Zyklu
s.
B
Umgebung und Applikationen
B.1 Entwicklungs- und Testumgebung
Die Implementierung und das Debugging der Applikationen und Scripte, sowie die
Evaluation der Ergebnisse erfolgte in folgender Umgebung:
• Ubuntu 10.04: Linux ubuntu 2.6.32-23-generic; 37-Ubuntu SMP Fri Jun 11 07:54:58
UTC 2010 i686 GNU/Linux
• Java: Version 1.6.0 18; OpenJDK Runtime Environment (IcedTea6 1.8) (6b18-
1.8-0ubuntu1); OpenJDK Client VM (build 14.0-b16, mixed mode, sharing)
• Python: Python 2.6.5 (r265:79063, Apr 16 2010, 13:09:56)
• Bash: 4.1.5(1)-release
• VMware Workstation: Version 7.1.0 build-261024
• Eclipse: Version 3.5.2, Build Id: M20100211-1343 (Galileo)
• HMModeler: Version 1.0.2
115
B. Umgebung und Applikationen 116
B.2 grepseq: Optionsbeschreibung
seqgrep Usage:
java -jar seqgrep.jar [OPTIONS] file
Options:
Mandatory arguments to long OPTIONS are mandatory for short
OPTIONS too.
-b REGEX Record start PATTERN
-c FIELD [REGEX] SCORETYPE [SCBOUND]
Assign the FIELD to select the score to
calculate the SCORETYPE
-e REGEX Record end PATTERN
-f FILETYPE Assign the TYPE of input file to parse
-h, --help Display this help and exit
-l Print order number with output lines
-m FIELD REGEX Assign the filter FIELD to match with REGEX
-n NUM Output the first/last NUM records
-o OUTFILE Place the output into OUTFILE
-q, --quiet, --silent Suppress all normal output
-r Reverse the result of sort comparisons
-s FIELD [REGEX] Assign the FIELD to sort
-t OP VALUE Assign the score threshold conditional
OPerator and VALUE
-x EXFILE Exclude records has lines start with strings
in EXFILE
-y, --contrary Invert the sense of matching
--version Output version information and exit
File type:
Specify the type of file to parse.
TYPE File Type Description
-------------------------------
F[ASTA] FASTA File Type
H[MMOO] HMMOO File Type (default)
↪→
B. Umgebung und Applikationen 117
Field control:
Specify field names to sort records and search for the pattern
in each record.
FIELD Line starts with... DataType
-------------------------------------------------------
_USR USER_PATTERN String
_CLC * Double
SQN "Sequenzbeschreibung: >" String
SQS "Sequenz:" String
PSC "Plain Score:" Double
PSF "Plain Score (Forward):" Double
SSC "Simple Score:" Double
SSF "Simple Score (Forward):" Double
RSC "Reverse Score:" Double
RSF "Reverse Score (Forward):" Double
SCS "Simple Corrected Score:" Double
SCSF "Simple Corrected Score (Forward):" Double
RCS "Reverse Corrected Score:" Double
RCSF "Reverse Corrected Score: (Forward):" Double
HCS "HMMer Corrected Score:" Double
HCSF "HMMer Corrected Score (Forward):" Double
SCORETYPE type:
Specify the type of score to calculate.
TYPE Score Type Description
------------------------------------------------------------------
NCDF NormCDF based on Mean (arithm. Average)
NCDFM NormCDF based on Median
NCDFT NormCDF based on Mean without the scores which are
out of bounds. CBOUND specifies the thresholds for
Least/Top scores.
REGular EXpression:
Specify regular expression to select the sort and search field.
For further help about JAVA compatible regular expressions visit:
http://java.sun.com/docs/books/tutorial/essential/regex/literals.html
Author
Roland J. Graf, University of Applied Sciences, Salzburg
Report bugs to: roland.graf@fh-salzburg.ac.at
Salzburg/Austria, August 2010
B. Umgebung und Applikationen 118
B.3 amodseq: Optionsbeschreibung
AmodSeq Usage:
java -jar AmodSeq.jar [OPTIONS] -p hmmfile -a alignedfile -o outfile
Options:
Mandatory arguments to long OPTIONS are mandatory for short OPTIONS too.
-a FILE Use the HMM model aligned sequences FILE
-f FILE Place the FASTA output into FILE
-d FILE Place further informations about the sequence
alignment process into FILE
-h, --help Display this help and exit
-l DIR Assign the location to search the external
sequence alignment application
-n NUM FILE Export first NUM chars of MSA sequence names
into FILE
-o FILE Place the output into FILE
-p FILE Use the HMM model parameter FILE created
with HMModeler
-t DIR Assign the location to create temporary files
-q, --quiet, --silent Suppress all normal output
-x NUM PARAMS Pass through NUM PARAMeters to the external
sequence alignment application
--version Output version information and exit
AUTHOR
Roland J. Graf, University of Applied Sciences, Salzburg
Report bugs to roland.graf@fh-salzburg.ac.at
Salzburg/Austria, July 2010
B. Umgebung und Applikationen 119
B.4 amodseq.ini: ClustalW Defaulteinstellungen
[XAPPARGS]
# specify path and application to improve MSA
/home/student/MT/Evaluation04-Roland/jarfiles/clustalw2
# specify the fixed parameters
# $$PROFILE1$$ = placeholder for fixed profile FASTA file
# $$PROFILE2$$ = placeholder for variable profile FASTA file
-ALIGN
-QUIET
-SEQUENCES
-OUTPUT=FASTA
-PROFILE1=$$PROFILE1$$
-PROFILE2=$$PROFILE2$$
-OUTFILE=$$OUTFILE$$
C
Daten- und Ergebnisdateien
C.1 Astral DB-Datei im FASTA Format
Das textbasierte FASTA-Format in ein in der Bioinformatik gangiges Dateiformat
zur Beschreibung der Primarstruktur von Proteinen. Eine Sequenzbeschreibung wird
mit einer Kopfzeile eingeleitet, die immer mit dem Zeichen ’>’ beginnt. Der In-
halt der Kopfzeile startet mit einem eindeutigen Namen oder einer ID. Die Ami-
nosauresequenzen selbst werden durch einen One-Letter -Code dargestellt. Kommentar-
zeilen werden mit dem Zeichen ’;’ eingeleitet.
Das nachfolgende Listing zeigt einen Auszug aus einer Astral-Datenbankdatei im FASTA
Format:
>d1dlwa_ a.1.1.1 (A:) Protozoan/bacterial hemoglobin {Ciliate (Para ↪→mecium caudatum) [TaxId: 5885]}
slfeqlggqaavqavtaqfyaniqadatvatffngidmpnqtnktaaflcaalggpnawt
grnlkevhanmgvsnaqfttvighlrsaltgagvaaalveqtvavaetvrgdvvtv
>d1s69a_ a.1.1.1 (A:) Protozoan/bacterial hemoglobin {Cyanobacteria ↪→(Synechocystis sp.), pcc 6803 [TaxId: 1143]}
stlyeklggttavdlavdkfyervlqddrikhffadvdmakqrahqkafltyafggtdky
dgrymreahkelvenhglngehfdavaedllatlkemgvpedliaevaavagapahkrdv
lnq
>d1idra_ a.1.1.1 (A:) Protozoan/bacterial hemoglobin {Mycobacterium ↪→tuberculosis, HbN [TaxId: 1773]}
gllsrlrkrepisiydkiggheaievvvedfyvrvladdqlsaffsgtnmsrlkgkqvef
faaalggpepytgapmkqvhqgrgitmhhfslvaghladaltaagvpsetiteilgviap
lavdvts
120
C. Daten- und Ergebnisdateien 121
C.2 HMMOO Alignment Result Datei
Das folgende Listing zeigt in gekurzter Form1 den grundsatzlichen Aufbau eines result-
Datensatz an, wie von der Applikation hmmoo generiert wird:
Sequenzbeschreibung: >d1dlwa_ a.1.1.1 (A:) Protozoan/bacterial hemo ↪→globin {Ciliate (Paramecium caudatum) [TaxId: 5885]}
Sequenz: slfeqlggqaavqavtaqfyaniqadatvatffngidmpnqtnktaaflcaalggpna ↪→wtslfeqlggqaavqavtaqfyaniqadatvatffngidmpnqtnktaaflcaalggpnawtgrnlk ↪→evhanmgvsnaqfttvighlrsaltgagvaaalveqtvavaetvrgdvvtv
Plain Score: -576.7240977002209
Plain Score (Forward): -569.0899435594409
Simple Score: -647.0919876918512
Simple Score (Forward): -647.0919876918512
Reversed Score: -581.9745888906357
Reversed Score (Forward): -572.879493576253
Simple Corrected Score: 70.36788999163025
Simple Corrected Score (Forward): 78.00204413241022
Reverse Corrected Score: 5.2504911904147775
Reverse Corrected Score (Forward): 3.789550016812086
HMMer Corrected Score: 65.11109477813987
HMMer Corrected Score (Forward): 71.77954943225613
BackTrack: [B, IB, IB, IB, IB, IB, IB, IB, IB, IB, IB, IB, IB, IB, ↪→...
IB, IB, IB, IB, IB, IB, IB, IB, IB, FB, M0, M1, M2, M3, M4, M5, M6, ↪→...
M84, M85, M86, M87, M88, M89, M90, M91, M92, M93, M94, M95, FE, IE, ↪→IE, IE, IE, IE, IE, IE, IE, IE, IE, E]
MSA:
SLFEQLGGQAAVQAVTAQFYANIQADATVATFFNGIDMPNQTNKTAAFLCAALGGPNAWTSLFEQLG ↪→GQAAVQAVTAQFYANIQADATVATFFNGIDMPNQTNKTAAFLCAALGGPNAWTGRNLKEVHANMGVS ↪→NAQFTTVIGHLRSALTGAGVAAALVEQTVAVAETVRGDVVTV
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIMM ↪→MMMMMMMMMMMMMMMMMMMMMMMMMMMIIMMMMMMMMMMMMMMMMMMMIMMIMMMMMMMMMIMMMMM ↪→MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMIIIIIIIIII
1Aus Platzgrunden wurde Teile des Datensatzes entfernt. Die fehlenden Zeilen wurden durch ein’...’ ersetzt.
C. Daten- und Ergebnisdateien 122
C.3 Z-Score Auswertung
Die folgenden Listings zeigen die vom Python-Tool evaluateResults.py generierte Aus-
wertung der Reverse Corrected Z-Scores fur die Zyklen 1 und 3. Bei den Z-Scores
handelt es sich um auf Other-Scores standardisierten Scores. Die Werte geben an, wie
groß die Distanz der Mittelwerte vom Other-Score zu den Mittelwerten der Family-
(Fam), Superfamily (SFam) und Fold-Scores (Fold) ist. Die Distanz wird in Vielfache
der Standardabweichung der Other-Scores angegeben. Die Medium Rank Werte geben
an, wie die Sequenzen im Durchschnitt gereiht wurden, wenn die Daten nach Scores
sortiert vorliegen. DBMember Werte zeigen die Vorkommen in der Testdatenbank.
C.3.1 Z-Score Auswertung Zyklus 1
------- Z-Score ------ --- Medium Rank ---- DBMembers
Family Fam SFam Fold Fam SFam Fold Fa -Sf -Fo
a.1.1.2 15.097 3.051 0.685 12.5 252.0 2580.5 26 11 4
a.138.1.1 16.786 2.913 nan 4.5 82.0 nan 10 12 0
a.25.1.2 17.574 0.531 -0.129 5.0 2828.0 5522.0 11 25 7
a.3.1.4 20.458 4.439 nan 2.5 23.0 nan 6 29 0
b.121.4.1 16.038 0.984 0.178 6.5 1232.0 3509.0 14 22 24
b.122.1.1 10.158 0.755 nan 3.5 2488.5 nan 8 20 0
b.45.1.1 8.812 -0.001 0.872 8.5 5335.5 660.0 18 6 3
b.47.1.2 39.318 4.695 nan 11.5 34.5 nan 24 20 0
b.6.1.1 15.484 1.959 -0.227 4.5 656.0 5485.0 10 35 2
c.31.1.3 20.659 1.929 nan 4.5 467.0 nan 10 10 0
c.36.1.9 31.904 2.457 nan 4.5 159.0 nan 10 20 0
c.45.1.2 47.279 5.085 nan 3.5 15.5 nan 8 12 0
c.46.1.2 11.496 5.516 nan 6.5 17.0 nan 14 7 0
c.72.1.1 30.190 2.473 1.119 6.0 171.5 943.0 13 8 9
d.131.1.2 9.155 0.915 nan 7.5 847.0 nan 16 6 0
d.14.1.5 19.021 -0.312 nan 3.5 6190.0 nan 8 26 0
d.153.1.4 24.103 -0.217 2.203 9.0 5431.0 195.0 19 17 1
d.32.1.2 19.039 3.099 nan 3.0 112.5 nan 7 24 0
d.38.1.5 12.584 3.298 nan 5.5 63.0 nan 12 33 0
e.3.1.1 17.904 5.282 nan 11.0 25.5 nan 23 4 0
Summary -------------------------------------------------------------
20.153 2.443 0.671
C. Daten- und Ergebnisdateien 123
C.3.2 Z-Score Auswertung Zyklus 3
------- Z-Score ------ --- Medium Rank ---- DBMembers
Family Fam SFam Fold Fam SFam Fold Fa -Sf -Fo
a.1.1.2 14.589 3.132 0.450 12.5 533.0 2727.5 26 11 4
a.138.1.1 10.422 0.983 nan 4.5 2290.0 nan 10 12 0
a.25.1.2 17.574 0.531 -0.129 5.0 2828.0 5522.0 11 25 7
a.3.1.4 18.059 3.814 nan 2.5 44.0 nan 6 29 0
b.121.4.1 16.633 1.097 0.178 6.5 1311.0 3787.0 14 22 24
b.122.1.1 8.941 0.757 nan 4.5 2629.0 nan 8 20 0
b.45.1.1 8.812 -0.001 0.872 8.5 5335.5 660.0 18 6 3
b.47.1.2 35.039 5.397 nan 11.5 36.5 nan 24 20 0
b.6.1.1 13.716 4.109 -0.009 5.5 133.0 4788.0 10 35 2
c.31.1.3 20.659 1.929 nan 4.5 467.0 nan 10 10 0
c.36.1.9 30.275 2.654 nan 4.5 112.5 nan 10 20 0
c.45.1.2 37.875 3.536 nan 3.5 53.5 nan 8 12 0
c.46.1.2 10.794 5.200 nan 6.5 18.0 nan 14 7 0
c.72.1.1 24.932 3.423 0.432 6.0 370.0 4388.0 13 8 9
d.131.1.2 9.155 0.915 nan 7.5 847.0 nan 16 6 0
d.14.1.5 19.021 -0.312 nan 3.5 6190.0 nan 8 26 0
d.153.1.4 23.352 -0.255 2.195 9.0 5161.0 210.0 19 17 1
d.32.1.2 17.504 2.964 nan 3.0 128.0 nan 7 24 0
d.38.1.5 12.202 4.187 nan 8.0 32.0 nan 12 33 0
e.3.1.1 17.016 10.016 nan 11.0 21.0 nan 23 4 0
Summary -------------------------------------------------------------
18.328 2.704 0.570
D
Quelltexte
D.1 amodseq: Alignment voralignierter Sequenzen
in das Profil-HMM
public int AlignSequenceToModel( ModelSequences modelSeq, ↪→AlignedSequence algSeq ) ↪→
throws AlignmentException
{
String seqID = algSeq.getName();
int seqOffset = 0;
int modelOffset = 0;
int matchOffset = 0;
int insertlen, endins;
StringBuilder seq = new StringBuilder(50000);
StringBuilder mask = new StringBuilder(50000);
seq.append(algSeq.getSequenceString());
mask.append(algSeq.getMaskString());
for(int i=0; i < mask.length(); i++)
{
switch( mask.charAt(i) )
{
case ’D’:
case ’M’:
matchOffset = modelSeq.getNextMatch(modelOffset);
if( matchOffset == modelOffset ){
modelOffset++;
seqOffset++;
} else {
insertlen = matchOffset - modelOffset;
mask.insert(i, StringUtils.RepeatChar(gapChar, insertlen));
124
D. Quelltexte 125
seq.insert(i, StringUtils.RepeatChar(gapChar, insertlen));
modelOffset = matchOffset+1;
i += insertlen;
}
break;
case ’I’:
for( endins=i+1; endins<mask.length()
&& mask.charAt(endins)==’I’; endins++);
matchOffset = modelSeq.getNextMatch(modelOffset);
if( matchOffset < 0 )
matchOffset = multipleAlignedSequences.getModelWidth();
insertlen = matchOffset - modelOffset;
if (insertlen > 0 ){
mask.insert(i, StringUtils.RepeatChar(gapChar, insertlen));
seq.insert(i, StringUtils.RepeatChar(gapChar, insertlen));
}
modelSeq.insertAtColumn(matchOffset,
StringUtils.RepeatChar(gapChar, endins - i).toString());
modelOffset = matchOffset + (endins - i);
i += insertlen + (endins - i - 1);
break;
default: // throws an exception because a wrong mask state
}
}
insertlen = modelSeq.getModelWidth() - seq.length();
if( insertlen > 0 )
{
seq.append(StringUtils.RepeatChar(gapChar, insertlen));
mask.append(StringUtils.RepeatChar(gapChar, insertlen));
}
algSeq.replace(seq.toString(), mask.toString());
return seq.length();
}
E
Datentrager
126