lab Grundlagen der Bioinformatik mod - Benutzer-Homepageskrausk/Bioinformatik/Bioinfo... · BLAST...

© modlab® 2007

modlab

Grundlagen der Bioinformatik

Ergänzende Unterlagen zu ausgewählten Kapiteln

Prof. Dr. G. Schneider

erstellt u.M.v. Markus Hartenfeller, Jan A. Hiß, Michael Meissner,

Ewgenij Proschak, Michael Schmuker, Yusuf Tanrikulu, Martin Weisel, Matthias Wirth

© modlab® 2007

Inhaltsverzeichnis

INHALTSVERZEICHNIS...................................................................................................... 3

SEQUENZALIGNMENT........................................................................................................ 5

PAARWEISES ALIGNMENT ....................................................................................................... 7 HEURISTIKEN FÜR LOKALES ALIGNMENT .............................................................................. 10 SIGNIFIKANZ LOKALER ALIGNMENTS .................................................................................... 14 MULTIPLES ALIGNMENT........................................................................................................ 15

Heuristiken für multiple Alignments ................................................................................ 17 LITERATUR ............................................................................................................................ 20

AUTOKORRELATIONSVEKTOREN............................................................................... 21

MATHEMATISCHE DEFINITION .............................................................................................. 22 TRANSFER IN DIE CHEMIEINFORMATIK.................................................................................. 23

Beispiel ............................................................................................................................. 25 ANWENDUNG VON KORRELATIONSVEKTOREN ...................................................................... 26 LITERATUR ............................................................................................................................ 27

NIPALS-ALGORITHMUS ................................................................................................... 29

PRINZIP UND SCHRITTE DES ALGORITHMUS .......................................................................... 30 LITERATUR ............................................................................................................................ 31

GRUNDLAGEN GRAPHENTHEORIE ............................................................................. 33

Graphen............................................................................................................................ 33 Subgraphen, Pfade, Zyklen, Cliquen................................................................................ 35

GRAPHENTHEORETISCHE PROBLEME UND IHRE LÖSUNGEN IN DER CHEMINFORMATIK......... 36 Subgraph-Isomorphismus................................................................................................. 36 Maximum Common Substructure ..................................................................................... 40 Cliquendetektion auf dem Assoziationsgraphen .............................................................. 41

LITERATUR ............................................................................................................................ 45

KÜNSTLICHE NEURONALE NETZE .............................................................................. 47

Genereller Aufbau von Feed-Forward-Netzen ................................................................ 47 DAS PERZEPTRON.................................................................................................................. 48

Das Single Layer Perzeptron ........................................................................................... 48 Die Lösung des XOR-Problems........................................................................................ 51

MEHRLAGIGE NETZE ............................................................................................................. 52 NETZTRAINING ...................................................................................................................... 53 PROBLEME UND LÖSUNGSANSÄTZE....................................................................................... 58 LITERATUR ............................................................................................................................ 61

PARTIKELSCHWARM-OPTIMIERUNG (PSO)............................................................. 63

DER ALGORITHMUS............................................................................................................... 64 ZUSAMMENFASSUNG ............................................................................................................. 66 LITERATUR ............................................................................................................................ 67

AMEISENALGORITHMEN................................................................................................ 69

AMEISENALGORITHMEN IN DER ANWENDUNG ...................................................................... 70 LITERATUR ............................................................................................................................ 72

Sequenzalignment 5

© modlab® 2007

Sequenzalignment

Die ersten Algorithmen zur Sequenzanalyse wurden in den 50er Jahren des 20. Jahrhunderts

erdacht, als die ersten Proteinsequenzen bekannt wurden. Mit der von Fred Sanger

entwickelten Methode zur Sequenzierung von DNA (Kettenabbruchreaktion, 1975) wurden

zudem schnell große Anzahlen an Nukleotidsequenzen verfügbar. Die hier vorgestellten

Alignment-Algorithmen dienen dem Vergleich und der Analyse von Sequenzalignments und

werden verwendet, um die strukturelle, funktionelle oder evolutionäre Verwandtschaft von

DNA- oder Proteinsequenzen aufzuklären.

Das Ausrichten zweier Sequenzen gelingt beim paarweisen Alignment durch

Gegenüberstellung von gleichartigen Positionen (Matches) oder möglichst ähnlichen

Positionen (Mismatches) und dem Einfügen von Lücken (Gaps) an den Stellen, an denen sich

die Sequenzen durch Punktmutationen (Substitutionen, Insertionen oder Deletionen)

unterscheiden. Dabei finden Scoring-Matrizen Anwendung, die überwiegend aus multiplen

Alignments bekannter homologer Sequenzen abgeleitet wurden und konservierten oder sehr

ähnlichen Positionen positive Werte (Scores) zuweisen, während unähnliche Positionen meist

negative Werte erhalten (es gibt auch Scoring-Matrizen, die auf physikochemischen

Eigenschaften beruhen und nicht aus Alignments abgeleitet wurden). Die Bestrafung von

Lücken erfolgt durch die Gap-Funktion γ (Formel 1), die negativ in den Score eines

Alignments eingeht. Dabei finden meist „affine“ Gap-Penalties Verwendung, bei denen ein

langes Gap weniger schlecht bewertet wird, als mehrere kurze. Dies gelingt dadurch, dass die

erste Stelle eines Gaps der Länge g mit einer Gap-Open-Penalty d stärker bestraft wird als die

restlichen Positionen der gleichen Lücke, die mit der Gap-Extension-Penalty (Formel 1) e

bestraft werden (d > e).

Sequenzalignment 6

© modlab® 2007

γ(g) = d + (g-1) · e. (1)

(Die lineare Gap-Penalty kann als Sonderfall der affinen Gap Penalty angesehen werden ,

wobei d = e.)

Bewertung eines Alignments am Beispiel:

Für die beiden DNA-Sequenzen TAAACAC und TGCATAC schlagen wir folgendes Alignment

vor:

TAA__ACAC

T__GCATAC

Wenn wir jeden Match mit +1 belohnen, jeden Mismatch mit –m und jeden Gap mit obiger

affiner Gap-Penalty γ bestrafen, erhält unser Alignment einen Score von 4 – m –2d – 2e. Wie

bereits erwähnt, werden üblicherweise Scoring-Matrizen verwendet, die Matches positive

Scores zuweisen, während Mismatches meist negative Werte erhalten. Dabei sorgt eine

Konvention dafür, dass ein Mismatch nicht strenger bestraft wird als zwei einzelne Gap-

Positionen, durch die man den Mismatch vermeiden könnte.

Während das optimale paarweise Alignment in Laufzeit O(n²) relativ schnell (aber oft immer

noch zu langsam für eine effiziente Datenbanksuche) berechnet werden kann, ist das optimale

multiple Alignment mit O(nk) (k = Anzahl der Sequenzen) nicht mehr effizient zu berechnen.

Hier finden Heuristiken (CLUSTAL W, Feng & Doolittle Ansatz) Verwendung, die

paarweises Alignment mit phylogenetischen Methoden kombinieren. So werden von den

Sequenzen ausgehend phylogenetische Bäume erstellt und die einander ähnlichsten

Sequenzen in paarweisen Alignments vereinigt, die wiederum zum fertigen multiplen

Alignment zusammengesetzt werden.

Im nun folgenden Abschnitt sollen Algorithmen zum paarweisen Alignment betrachtet

werden.

Sequenzalignment 7

© modlab® 2007

Paarweises Alignment

Globales Alignment

Wir haben zwei Sequenzen x,y der Länge n,m gegeben und vermuten, dass diese miteinander

verwandt sind. Um nicht alle exponentiell viele Alignments berechnen zu müssen, benutzen

wir einen dynamischen Programmierungsansatz: Der Needleman-Wunsch-Algorithmus

(1970) liefert uns optimale globale Alignment in Laufzeit O(n · m). Dazu erstellen wir uns

einen (n × m) Alignmentgraphen, in dem wir in allen Feldern je drei Kanten beschriften (Abb.

1).

Die mit Di,j beschrifteten diagonalen Kanten repräsentieren einen Match bzw. einen Mismatch

zwischen den jeweiligen Positionen i,j in den Sequenzen x,y und werden mit den zugehörigen

Scores S beschriftet. Die vertikalen bzw. horizontalen Kanten (V, H) repräsentieren jeweils

einen Gap in der ersten bzw. zweiten Sequenz. Indem wir uns bereits errechnete Scores der

besten Sequenzanfänge merken und mit dynamischen Programmieren weitergeben, erhalten

wir die o.g. Laufzeitbeschleunigung. So trägt jede mit Di,j beschriftete Kante den Score des

besten Alignments bis zu dieser Stelle (Formel 2, hier ausgehend von Di+1,j+1).

Di+1,j+1 = S(xi+1, yj+1) + max {Di,j + Vi,j + Hi,j} (2)

x1

y1

x2 xn

y2

ym

D1,1

H1,0

V0,1

Hn,m

Dn,m Vn,m

Di,j Vi,j

Hi,j

xi xi+1

yj

yj+1 Di+1,j+1Vi,j+1

Hi+1,j

Abbildung 1: Die Kanten des Alignmentgraphen tragen den Score des besten Alignments bis zu dieser Stelle.

Sequenzalignment 8

© modlab® 2007

Bei den Rekursionen für horizontale und vertikale Kanten werden keine Scores verteilt, da

Gaps keine Beiträge zum Alignment liefern. An dieser Stelle fließen jedoch die Gap-Open-

oder Gap-Extension-Penalties (d,e) ein (Formel 3).

Hi+1,j = max {Di,j - d, Vi,j - d, Hi,j - e} (3)

Nachdem wir alle Kanten beschriftet haben, erhalten wir den Score des optimalen globalen

Alignments (Formel 4).

S* = max {Dn,m, Vn,m, Hn,m} (4)

Den dazugehörigen Pfad erhalten wir per Traceback über die maximalen Kanten. Das

bedeutet wir laufen die Kanten rückwärts bis zum Start und entscheiden uns jeweils für die

Kante, die den maximalen Beitrag zum Score geliefert hat (Abb. 2).

Abbildung 2: Der Traceback über die Kanten mit den maximalen Beschriftungen liefert den Pfad zum optimalen Alignment.

Ein Wort zur Laufzeit: Wir beschriften insgesamt 3·n·m Kanten und benötigen für den

Traceback maximal n+m Kanten und bleiben mit 3·n·m +n+m asymptotisch in O(n·m).

Lokales Alignment

Wenn wir vermuten, dass nur ein Teil (xk,...,xi) einer Sequenz x mit einem Teil (yl,...,yj) einer

Sequenz y übereinstimmt, verwenden wir den Smith-Waterman-Algorithmus (1981), der das

optimale lokale Alignment für die beiden Sequenzen in O(n·m) findet. Diese Suchstrategie

wird z.B. verwendet, wenn wir eine Sequenz gegen eine ganze Datenbank alignieren möchten

oder überprüfen wollen, ob z.B. ein Hefeprotein irgendwo im Humangenom vorkommt.

S*

Sequenzalignment 9

© modlab® 2007

Wir benutzen wieder die affine Gap-Penalty und beschriften die Kanten Di,j, Hi,j, Vi,j wie beim

Needleman-Wunsch-Algorithmus, jedoch mit folgendem Unterschied: Wir geben jeder Kante

die Möglichkeit, ein neues lokales Alignment mit Score ≥ 0 zu starten. Das bedeutet, dass

eine Kante nicht ein möglicherweise negatives Erbe von schlechten Sequenzanfängen

weiterführen muss, sondern stets ein neues lokales Alignment starten kann. Mit anderen

Worten: Die Teile die außerhalb des optimalen lokalen Alignments liegen, gehen nicht

negativ in den Score ein (Abb. 3).

Abbildung 3: Das optimale lokale Alignment – kein anderes lokales Alignment erreicht einen größeren Score.

Die Rekursionen ändern sich wie folgt (Formeln 5, 6, 7):

Di+1,j+1 = max {0, Di,j + S(xi+1, yj+1), Hi,j + S(xi+1, yj+1), Vi,j + S(xi+1, yj+1)} (5)

Hi+1,j = max {0, Di,j - d, Hi,j - e, Vi,j - d} (6)

Vi,j+1 = max {0, Di,j - d, Hi,j - d, Vi,j - e} (7)

Das optimale lokale Alignment erhalten wir per Traceback über die Kanten, die den

maximalen Beitrag geliefert hatten. Dabei haben wir uns während des Beschriftens der

Kanten den maximalen Score gemerkt und gehen von dieser Kante (nicht von Dn,m, Hn,m oder

Vn,m wie beim Needleman-Wunsch-Algorithmus) nach links oben vor, bis wir eine Kante

erreichen, die den Score 0 trägt. Der Teil zwischen dem Score 0 und dem maximalen Score ist

unser lokal optimales Alignment.

Laufzeit und Speicherplatz betragen wie beim Needleman-Wunsch-Algorithmus O(n·m), was

deutlich besser ist, als exponentielle Laufzeit, aber für die Suche in großen Datenbanken zu

xk... ...xi

yl

yj

...

...

Sequenzalignment 10

© modlab® 2007

langsam ist. Hier finden die Heuristiken BLAST oder FASTA Anwendung, die zwar nicht

unbedingt das optimale lokale Alignment finden, dafür aber deutlich schneller laufen und

meist gute Ergebnisse liefern.

Heuristiken für lokales Alignment

BLAST (Basic Local Alignment Search Tool)

BLAST ist eine Heuristik für lokale Alignments, die mit einer Laufzeit von O(n+m) deutlich

schneller ist als der Needleman-Wunsch-Algorithmus, jedoch nicht in jedem Fall das optimale

lokale Alignment findet.

Gegeben seien eine kürzere Sequenz S der Länge m und eine lange Sequenz T mit Länge n,

gegen die S in einem (oder mehreren) lokalen Alignment(s) aligniert werden soll. Dazu wird

S in kleine Teilwörter zerteilt (Default: Wortlänge |w| = 3 bei Proteinen, |w| = 11 bei

Nukleinsäuren). Hierbei fügt BLAST auch ähnliche Wörter hinzu, deren Score einen

bestimmten Schwellenwert t überschreiten (diese leicht veränderten Suchwörter stehen für

Mutationen, die sich während des Evolutionsprozesses eingeschlichen haben könnten). Aus

diesen Teilwörtern lässt sich ein Schlüsselwortbaum so konstruieren, dass die lange Sequenz

T fast in linearer Zeit durch den Suchbaum aus den Fragmenten des Suchstrings S geschoben

werden kann. Dies gelingt dadurch, dass die Suchmuster so in den Baum eingefügt werden,

dass sie sich erst aufspalten, wenn sie sich unterscheiden. Z.B. für Sequenz S = KLLGV

(Abb. 4).

Abbildung 4: Schlüsselwortbaum für Sequenz S.

L VG

L L G

K L

Sequenzalignment 11

© modlab® 2007

So kann man im Fall eines Mismatches bei der Mustersuche direkt an die richtige Stelle eines

anderen Musters mit gleichem Präfix springen. Dies bringt die gewünschte

Laufzeiteinsparung (Vergleiche hierzu auch weiterführende Literatur. Stichwort: „Aho-

Corasick-Algorithmus“).

Wir schieben nun den langen Text T durch den Suchbaum und geben alle Matches aus (alle

Fragmente mit Score ≥ t). Diese lassen sich in einem Alignmentgraphen (Dotplot) darstellen

und repräsentieren Stücke, die in beiden Sequenzen identisch (bzw. ähnlich im Falle der

zusätzlich eingefügten Strings mit Score ≥ t) sind. In einem weiteren Schritt verlängern wir

diese gefundenen Alignments entlang der Diagonalen zu sogenannten HSPs (High Scoring

Segment Pairs). Diese sind so definiert, dass sie sich durch Verlängern oder Verkürzen in

ihrem Score nicht mehr verbessern lassen. Wir können uns nun eine Reihe guter lokaler

Alignments ausgeben lassen, die einen bestimmten Mindestscore überschreiten oder einfach

das MSP (Maximum Scoring Segment Pair) betrachten. Wie schon angedeutet findet BLAST

nicht unbedingt das optimale lokale Alignment, liefert aber meist sehr gute Ergebnisse. Die

Laufzeit ist mit O(n+m) deutlich schneller als beim Smith-Waterman-Algorithmus.

Anmerkung: In der klassischen BLAST Version werden HSPs nur entlang der Diagonalen

verlängert. Beim sog. GAPPED BLAST werden nur solche HSPs betrachtet, die (innerhalb

eines begrenzten Abstands zueinander) auf einer gemeinsamen Diagonalen liegen, dass heißt

einem gemeinsamen gaplosen Alignment angehören. Zwischen diesen HSPs suchen wir nun

nach örtlichen Exkursionen, die unsere beiden HSPs mit einem weiteren HSP über Gaps

verbinden kann (Abb. 5).

Abbildung 5: Gapped BLAST: einzelne HSPs werden über Gaps verbunden, wenn sich dadurch der Score verbessern lässt.

Sequenz 1

Sequenzalignment 12

© modlab® 2007

Um weiter entfernte Sequenzen zu finden empfiehlt sich PSI-BLAST (Position Sensitive

Iterated-BLAST). Hierbei wird aus signifikanten Treffern einer BLAST-Suche ein multiples

Alignment erstellt, aus dem ein Profil abgeleitet werden kann. Im Profil sind die jeweiligen

Wahrscheinlichkeiten für jedes Nukleotid bzw. jede Aminosäure für die jeweilige Position im

Alignment vermerkt. Dieses Profil wird nun zur Datenbanksuche eingesetzt, die

resultierenden signifikanten Treffer wieder aligniert und ein neues (noch variableres) Profil

erstellt. Nach mehreren Iterationen endet das Verfahren. PSI-BLAST findet auch weiter

entfernte Sequenzen (auch < 30% Identität).

Die wichtigsten Datenbanken (Genbank, EMBL und DDBJ) lassen sich allesamt mit BLAST

durchsuchen. Die zahlreichen BLAST-Varianten sind in der BLAST-Suite zusammengefasst:

http://www.ncbi.nlm.nih.gov/blast/producttable.shtml#pstab

Die wichtigsten Varianten:

blastn: zur Suche in Nukleotidsequenzdatenbanken.

blastp: zur Suche in Proteinsequenzdatenbanken.

blastx: Vergleicht eine Nukleotidsequenz (in allen Leserastern

translatiert) gegen eine Proteindatenbank.

tblastn: Vergleicht eine Proteinsequenz gegen eine Nukleotiddaten-

bank (dynamisch in allen Leserastern translatiert)

tblastx: Vergleicht die six-frame-Translation einer Nukleotidsequenz

gegen die six-frame-Translationen einer Nukleotidsequenz-

datenbank.

FASTA

Auch hier soll ein kurzer Text S in einem langen Text T gefunden werden. Wir gehen ähnlich

wie bei BLAST vor und zerteilen Text S in kleinere Bruchstücke. Der Parameter „ktup“

bestimmt hier die Länge der Bruchstücke, die gematcht werden sollen (Default: ktup = 6 für

DNA, 2 für Proteine). Da die Bruchstücke hier kleiner sind als bei BLAST, suchen wir nicht

mit Pattern-Matching-Algorithmen nach Übereinstimmungen, sondern mit dem sog. 'lookup

Sequenzalignment 13

© modlab® 2007

table', einer Hashing-Tabelle. Im Lookup-Table wird das Auftreten aller Aminosäuren bzw.

Nukleotide in der Sequenz T aufgeführt, z.B. für die Sequenz T= A1G2T3A4C5A6 (Tab. 1).

Tabelle 1: Lookup-Table für Sequenz T

Beginnt ein Suchstring der Länge ktup aus S mit einem 'A', erkennen wir im lookup-table,

dass im Text T ein 'A' an den Positionen 1,4 und 6 auftritt. Steht in S an zweiter Postion ein

'G', prüfen wir, ob es ein 'G' in T an den Positionen 2,5 oder 7 gibt (also ein Zeichen nach den

schon gefundenen 'A's.), usw. Man kann so sehr schnell alle "hotspots" aufspüren. Das sind

genau die Stellen, an denen T und S an mindestens ktup Positionen übereinstimmen. Die

Übereinstimmungen tragen wir in einen Dotplot ein.

Wir können nun die Init1-Scores berechnen, d.h. wir verlängern alle hotspots zu HSPs (s.o.)

entlang der Diagonalen. Danach berechnen wir die Initn-Scores, dass heißt, wir versuchen

lokal die besten n HSPs durch Gaps zu verbinden, s.d. ein noch höherer Alignment-Score

resultiert. Die Breite der Nachbarschaftssuche wird dabei auf einen bestimmten Rahmen

beschränkt: Ausgehend von den Init1's betrachten wir die umliegenden 32 Diagonalen (16

höher, 16 tiefer).

Ein besonderes Merkmal von FASTA ist die anschließende "Exakte Scorebestimmung": Für

den Bereich in dem das Alignment mit dem höchsten Initn-Score wird nochmal ein exaktes

lokales Alignment mit dem Smith-Waterman-Algorithmus berechnet. Dabei ist der

betrachtete Rahmen noch etwas größer als der Bereich in dem das beste Initn-Alignment liegt,

um eventuelle Fehler von FASTA auszugleichen (Abb. 6).

Fazit: BLAST vs. FASTA

BLAST läuft schneller, weil die Suchkriterien hier strenger sind als bei FASTA (FASTA

findet mehr Treffer, wg. Default ktup = 2, während |w| = 3 bei BLAST (für Proteine)). Dafür

findet FASTA aber mit höherer Wahrscheinlichkeit das optimale Alignment, durch

Anwendung des Smith-Waterman-Algorithmus.

Symbol Position A 1,4,6 G 2 T 3 C 5

Sequenzalignment 14

© modlab® 2007

Suchraum um Init1's

Smith-Waterman

Sequenz 1

Abbildung 6: Exakte Scorebestimmung von FASTA.

Signifikanz lokaler Alignments

Die Signifikanz eines Alignments drückt aus, wie aussagekräftig dessen Score tatsächlich ist:

Ist der Score das Ergebnis eines zufälligen Prozesses oder ist der Score so hoch, dass wir

ausschließen können, dass die Ähnlichkeit der Sequenzen zufällig ist? In diesem Fall gehen

wir davon aus, dass die Sequenzen homolog sind. Hier sollen mit dem Z-Wert, P-Wert und

dem E-Wert drei Maße für Signifikanz vorgestellt werden.

Z-Wert

Ist der Score S unseres Alignments besser als der einer zufällig (aus den gleichen Bausteinen)

zusammengewürfelten Sequenz? Dazu berechnen wir üblicherweise 100 neue Alignments,

dazu den Mittelwert ihrer Scores und die Standardabweichung (Formel 8).

Z-Score von S = (S – Mittelwert) / Standardabweichung (8)

Ein Z-Wert von ≥ 5 gilt als signifikant. Dies tritt ein, wenn der beobachtete Score deutlich

höher ist, als der Mittelwert der zufällig erstellten Sequenzen. Je höher ein Z-Score, desto

geringer ist dabei die Wahrscheinlichkeit, dass der beobachtete Alignmentscore nur Zufall ist.

Ein Z-Score von 0 bedeutet, dass der beobachtete Score nicht größer ist als der einer

„zufälligen Sequenz“.

Sequenzalignment 15

© modlab® 2007

P-Wert

Hier berechnen wir die Wahrscheinlichkeit für das Ereignis, dass ein so hoher Score wie der

beobachtete Score S zufällig zustande kam (Nullhypothese). Fällt die Wahrscheinlichkeit für

alle Alignments der gleichen Länge, die einen Score S erreichen oder übertreffen unter das

Signifikanzniveau von 5%, verwerfen wir die Nullhypothese und gehen von Homologie der

Sequenzen aus. Das heißt, die Ähnlichkeit der betrachteten Sequenzen ergibt sich durch

gemeinsame Abstammung und nicht durch Zufall.

E-Wert

Der E-Wert (Formel 9) gibt den Erwartungswert für ein unabhängiges Alignment an, dass

dieses einen Score ≥ S erreicht.

E-Wert für S = n · m · k · e-λS . (9)

λ und k sind hier Skalierungsfaktoren für die Größe des Suchraums und die gewählte

Scoringmatrix.

Der E-Wert wird besonders häufig zur Bestimmung der Signifikanz von homologen

Alignments verwendet, da er (anders als Z- und P-Wert) Sequenzlängen berücksichtigt. Das

bedeutet, liefern zwei Sequenzen A,B den gleichen Score S und A ist kürzer als B, so hat B

einen deutlich höheren E-Wert (ist also weniger signifikant), da es ja mehr Positionen

benötigt hat, um den Score S zu erreichen. Die unterschiedlichen Sequenzlängen gehen mit

(n·m) in die Formel des E-Wertes ein.

Multiples Alignment

Das multiple Alignment liefert im Gegensatz zum paarweisen Alignment genauere

Informationen über Aminosäureverteilungen an bestimmten Positionen, da sich hieraus

ableiten lässt, wie stark konserviert eine Position tatsächlich ist. Darüberhinaus geben sie

Aufschluss darüber, wie wahrscheinlich der Austausch einer bestimmten Aminosäure durch

Sequenzalignment 16

© modlab® 2007

eine andere ist. Hierauf beruhen z.B. die PAM-Scoring-Matrizen, die sich aus den

Austauschwahrscheinlichkeiten von Aminosäuren aus Alignmentblöcken errechnen lassen.

Zudem lassen sich aus multiplen Alignments von Proteinen einer Familie Profile ableiten, mit

denen man nach weiteren Vertretern dieser Familie suchen kann.

Während für paarweises Alignment effiziente Algorithmen vorhanden sind, ist die

Berechnung des optimalen Alignments für mehr als zwei Sequenzen nicht mehr effizient

möglich. Hier müssten wir den Weg mit optimalem Score durch eine n-dimensionale Matrix

finden, was exponentielle Laufzeit benötigt. Abhilfe schaffen hier Heuristiken für multiple

Alignments.

Scores für multiples Alignment

Der Gesamtscore eines multiplen Alignments m errechnet sich als die Summe der Scores

seiner Spalten i. Zusätzlich gehen hier die Gap-Penalties G negativ ein (Formel 10).

∑i

i )S(m+G=S(m) (10)

Die Scores der einzelnen Spalten lassen sich auf verschiedene Weise berechnen:

Sum of Pair Scores

Hierzu berechnen wir die Scores aller n! möglichen Paare einer Alignmentposition (z.B. mit

PAM- oder BLOSUM-Matrizen). Ein Problem beim Sum of Pair Score ist, dass kleine

evolutionäre Ereignisse überbewertet werden. Ist nämlich nur eine Position einer sonst

konservierten Spalte mutiert, so geht der (meist negative) Score für den Mismatch mit allen

anderen (n-1) Positionen n-1 mal in den Spaltenscore ein.

Entropie-Scores

• Shannon-Entropie (Formel 11)

∑−|A|

=kkk pp=H

12log , wobei (11)

Sequenzalignment 17

© modlab® 2007

|A| = Die Länge des Alphabets (z.B. 20 Aminosäuren)

pk = Die relative Häufigkeit des Elementes k in der betrachteten Spalte

Die Shannon-Entropie ist ein Minimum-Entropie-Score, d. h. eine durchweg

konservierte Spalte (Optimum) erhält den Score 0, während der Score bei immer

größer werdender Unordnung (Entropie) immer größere (positive!) Werte annimmt

(das Maximum ist abhängig von der Größe des Alphabets).

Der sog. Information-Content I (Formel 12) kann als die „reverse“ Shannon-Entropie

aufgefasst werden. Er drückt aus, wie aussagekräftig eine Spalte im multiplen

Alignment ist (je höher der Wert, umso besser).

∑−|A|

=kkkact pp+|A|=HH=I

122max loglog . (12)

• Kullbach-Leibler-Distanz (Formel 13)

∑k k

kk

q

pp=Q)||H(P log (13)

Die Kullbach-Leibler-Distanz beschreibt die „Distanz“ zwischen der beobachteten

Häufigkeit pk (also der relativen Häufigkeit einer Aminosäure oder eines Nukleotids k

in der Spalte) und der Hintergrundverteilung qk (der erwarteten Wahrscheinlichkeit für

das Auftreten einer Aminosäure oder eines Nukleotids). Mit anderen Worten: Die

durch einen Mismatch in der Spalte hervorgerufene Entropie wird noch erhöht, wenn

die betreffende Position (AS oder Nukleotid) besonders selten auftritt.

Heuristiken für multiple Alignments

Die Berechnung des optimalen multiplen Alignments von n Sequenzen ist zeitaufwendig. So

ließe sich zwar in einem n-dimensionalen Ansatz ähnlich zum Needleman-Wunsch-

Algorithmus durch die Berechnung der jeweils besten Scores in einer n-dimensionalen Matrix

durch dynamisches Programmieren das optimale multiple Alignment berechnen, jedoch wäre

Sequenzalignment 18

© modlab® 2007

die Laufzeit immer noch exponentiell. Das Problem ist, das jede der n Positionen in einer

Spalte entweder ein Gap oder eine tatsächliche Position (AS oder Nukleotid) sein könnte.

Wenn wir verbieten, dass alle Positionen einer Spalte Gaps sein können, beschert uns dies

immer noch 2n-1 mögliche Kombinationen, deren Scores es über m Spalten hinweg zu

maximieren gilt. Es gibt Branch & Bound-Ansätze, die nicht die Sequenzanfänge aller

Alignments berechnen (wie im Brute-Force-Ansatz), sondern unwahrscheinliche Alignments

von der Berechnung ausschließen. Diese Methoden benötigen jedoch immer noch

exponentielle Laufzeit, auch wenn sie um eine Konstante besser sind. Nachfolgend sollen

zwei Heuristiken vorgestellt werden, die zwar einen Verlust der Genauigkeit in Kauf nehmen,

dabei aber deutlich schneller sind als Brute-Force-Ansätze und dennoch gute Ergebnisse

bringen.

Feng-Doolittle-Algorithmus

Vorgehen:

1. Aligniere alle Sequenzen paarweise mit dem Needleman-Wunsch-Algorithmus.

2. Weise allen n-1 Sequenzpaaren Distanzen zu. Der Feng-Doolittle-Algorithmus

berechnet hierzu den sog. „difference score“. Dieser entspricht in etwa dem Reversal

des erhaltenen Scores eines paarweisen Alignments, d.h. je höher der Score des

Alignments, desto niedriger der difference score. So wird eine hohe Ähnlichkeit

zwischen zwei Sequenzen in eine niedrige Distanz übersetzt.

3. Berechne aus den Distanzen einen phylogenetischen Baum mit UPGMA.

4. Aligniere längs der Kanten (ausgehend von den Blättern) Sequenz-Paare, Sequenz-

Block-Paare und Block-Paare (Abb. 7).

Seq 1 Seq 2 Seq 3 Seq 4 Seq 5

Block (1,2)

Block (4,5)

Block (1,2,3)

Block (1,2,3,4,5)

Abbildung 7: Multiples Alignment mit Feng-Doolittle: Aligniere zunächst Blätter untereinander, dann Blätter gegen Blöcke und schließlich Blöcke untereinander.

Sequenzalignment 19

© modlab® 2007

Zunächst werden je zwei Sequenzen zu einem Block zusammengefügt, und gegen diesen

Block können später weiter Sequenzen oder andere Blöcke aligniert werden. Ein einmal in

einem Block fixiertes Alignment bleibt für das multiple Alignment mit all seinen Gaps nach

der Regel „once gap, always gap“ erhalten. Mit einer Ausnahme: In den fixierten Block kann

ein globales Gap eingefügt werden, um eine weitere Sequenz gegen den Block zu alignieren.

Eine Sequenz, die neu zu einem Block hinzugefügt werden soll, wird zunächst mit allen

Sequenzen des Blocks paarweise aligniert (n-1 viele paarweise Alignments). Die Sequenz

wird dann gemäß des besten paarweisen Alignments in den Block eingefügt. Als Ausgabe

erhalten wir einen Block, der das Alignment aller unserer Sequenzen beinhaltet, also das

multiple Alignment.

CLUSTAL W

Clustal W verfolgt eine ähnliche Strategie wie Feng-Doolittle. Deshalb sollen hier nur die

wesentlichen Unterschiede vorgestellt werden. Zunächst werden auch hier alle paarweisen

Alignments berechnet und Distanzen zugewiesen. Dabei kennt CLUSTAL W neben dem

difference score aus dem Alignment Score wahlweise noch einen heuristischen Ansatz zur

Berechung der Distanzen. So kann die Distanz in einem beschleunigten Verfahren über die

Anzahl der identischen Treffer bei einer ktup-Suche (siehe FASTA) berechnet werden. Aus

den Distanzen wird wieder ein phylogenetischer Baum aufgebaut, hier aber nach der

Neighbor-Joining-Methode. Ausgehend von den Blättern werden nun wieder zunächst

Sequenzen miteinander aligniert. Diese werden aber nicht zu fixierten Blöcken

zusammengefasst, sondern an den Knoten zu Profilen geclustert. Ein Profil ist eine Matrix,

die für jede der 20 Aminosäuren (oder 4 Nukleotide) die relative Häufigkeit ihres Auftretens

an dieser Stelle im multiplen Alignment beinhaltet. Um zu verhindern, dass eng verwandte

Sequenzen das Profil in ihre Richtung verzerren, wird ihr Einfluss auf das Profil rechnerisch

herabgesetzt. Als Ausgabe erhalten wir ein Profil, aus dem sich ablesen lässt, welche Stellen

im Alignment hoch konserviert sind (evtl. funktionstragend) und welche eher nicht.

Sequenzalignment 20

© modlab® 2007

Literatur

Needleman SB, Wunsch CD (1970). A general method applicable to the search for similarities in the amino acid sequence of two proteins. J Mol Biol. 48(3): 443-53. Smith TF, Waterman MS (1981). Identification of common molecular subsequences. J Mol Biol. 147(1): 195-7. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ (1990). Basic local alignment search tool. J Mol Biol. 215(3): 403-10. Pearson WR (1994). Using the FASTA program to search protein and DNA sequence databases. Methods Mol Biol. 24: 307-31. Lesk A. (2002). Introduction to Bioinformatics. Oxford University Press, Oxford. Feng DF, Doolittle RF (1990). Progressive alignment and phylogenetic tree construction of protein sequences. Methods Enzymol. 183: 375-87. Thompson JD, Higgins DG, Gibson TJ (1994). CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22(22): 4673-80.

Autokorrelationsvektoren 21

© modlab® 2007

Autokorrelationsvektoren

Molekulare Deskriptoren sind im Wirkstoffdesign und für QSAR Modelle unentbehrlich,

um physikochemische Eigenschaften von Molekülen mit mathematischen Methoden zu

beschreiben. Sie ermöglichen es, umfangreiche Datenbanken mit Millionen von

Verbindungen digital abzuspeichern und chemoinformatische Methoden anzuwenden.

Autokorrelationsvektoren stellen nur eine Klasse der molekularen Deskriptoren dar. Durch die

Zusammenfassung von Korrelationswerten zu Vektoren lässt sich ein molekularer Deskriptor

erstellen, der Ähnlichkeitsberechnungen von Substanzen auf Basis der Verteilung von

Atomeigenschaften und deren Abständen zueinander zulässt. Hierbei erfolgt eine Abstraktion

der chemischen Bindungen zwischen den Atomen eines Moleküls und der Molekülstruktur im

kartesischen Raum.

Grundlegend gibt es zwei Anforderungen an einen Deskriptor. Zum einen soll er den

Vergleich von Molekülen zulassen, und zum anderen soll er praktisch anwendbar sein und

eindeutig aus den bekannten Daten berechnet werden können. Beides trifft auf

Korrelationsvektor-basierte Deskriptoren zu:

a) die Darstellung als Vektor ermöglicht eine einfache mathematische Quantifizierung

der Ähnlichkeit von Molekülen (durch Vektoralgebra), und

b) die minimale Anforderung zur Erstellung eines Korrelationsvektors ist die

zweidimensionale Molekülstruktur, die meist vorhanden ist oder leicht bestimmt

werden kann.

Bei der Verwendung von Autokorrelationsvektoren als molekularen Deskriptor sollten auch

die Nachteile stets bedacht werden. Erstens ist hier die Vernachlässigung der vollständigen

Molekülstruktur inklusive wichtiger dreidimensionaler Informationen zu nennen. Außerdem


© modlab® 2007

ist keine Rekonstruktion der Molekülstruktur anhand des zugehörigen Korrelationsvektors

möglich. Der Abstraktionsgrad des Deskriptors ist zu groß, als dass eine bijektive Darstellung

möglich wäre.

Mathematische Definition

Die Wurzeln der Autokorrelationsfunktionen liegen in der Funktionsanalyse von

physikalischen Methoden. Im Allgemeinen hat sie die Form:

∫ +⋅=b

a

t dxtxfxfAC )()( . (1)

Dabei stellt die Funktion f(x) eine beliebige Funktion dar, die auf dem Intervall [a, b+t]

definiert ist. Die Berechnung des Korrelationswertes an der Stelle t lässt sich folgendermaßen

beschreiben (Abbildung 1):

Gegeben seien zwei Punkte in dem Intervall [a, b+t], deren Abstand zueinander t beträgt.

Beide Punkte laufen nun unter Beibehaltung des Abstands über das Intervall. Für jeden

Zeitpunkt berechnet die Autokorrelationsfunktion das Produkt der Funktionswerte an beiden

Punkten und bildet das Integral über das gesamte Intervall. Als Ergebnis geht aus der

Rechnung ein Wert hervor, der nicht von einem Funktionswert an einer bestimmten Stelle x

abhängt, sondern vielmehr von der Verteilung der Funktionswerte im Intervall zu einem

gegebenen Zeitpunkt.

Abbildung 1: Darstellung von zwei möglichen Punkten im Abstand t einer gegebenen Funktion f(x), die autokorreliert werden können.


© modlab® 2007

Diese Tatsache ermöglicht vor allem die Anwendung von Autokorrelationsfunktionen in der

Analyse von physikalischen Spektren. Beispiele hierfür sind zum Beispiel NMR-

Spektrenauswertungen und Spracherkennungsprozesse. In den meisten Fällen ist zu einem

gegebenen Spektrum keine Funktion vorhanden, sondern es liegt nur der Graph vor. Die

Stärke der Autokorrelationsfunktion ist die Erkennung von sich periodisch wiederholenden

Signalen. Da die Produkte aus den Spitzen mit einem Abstand t in das Integral mit einfließen,

machen sich wiederkehrende Elemente mit einer Periodenlänge t durch Maxima in der

Autokorrelationsfunktion bemerkbar. Bei einer unterschiedlichen Periode wird der

Funktionswert an der Spitze der Funktion mit einem möglicherweise kleineren Wert

multipliziert, woraufhin die Wertigkeit des gesamten Integrals kleiner wird. Diese Eigenschaft

ist es, die die Verwendung von Autokorrelationsfunktionen insbesondere für das Filtern von

stark verrauschten Signalen auszeichnet.

Transfer in die Chemieinformatik

Moreau und Broto waren 1980 die ersten, die die Autokorrelationsfunktionen in die Welt der

Chemieinformatik einführten. Die Hauptintention war es, Aussagen über die Verteilung von

Atomeigenschaften, die als numerische Werte vorlagen, im Molekülgraphen machen zu

können. Der Transfer der klassischen Autokorrelations-funktionen auf ein als Strukturformel

vorliegendes Molekül wurde in folgender Weise vollzogen:

• Anstatt der Stelle x auf einem Funktionsgraphen wurden die Atome im Molekül als zu

betrachtende Punkte verwendet. Die Wasserstoffatome wurden hierbei nicht in

Betracht gezogen. Alle anderen Atome wurden beliebig durchnummeriert, um eine

eindeutige Identifizierung der Punkte zu ermöglichen. Das Intervall [a, b+t] wurde

also durch das gesamte Molekül dargestellt. Aufgrund der Tatsache, dass es nur eine

endliche Anzahl von Punkten (Atome im Molekül) gibt, und keine kontinuierliche

Funktion vorliegt, konnte das Integral der Funktion durch die Aufsummierung der

Korrelationswerte ersetzt werden.

• Als Funktionswerte f(x) wurden die Atomeigenschaften des Atoms x betrachtet. Meist

wurde die Funktion als p(x), wie property, bezeichnet. In jeder

Autokorrelationsfunktion wirdnur eine Atomeigenschaft mit sich selbst korreliert.


© modlab® 2007

• Schließlich wurde der Abstand t als die minimale Anzahl der Bindungen zwischen

zwei Atomen definiert (topologische Distanz). Hierbei wird nicht zwischen Doppel-

und Einfachbindungen unterschieden.

•

Die auf diese Weise modifizierte Autokorrelationsfunktion entspricht der Formel:

∑∈

⋅=)(,

)()(dMji

d jpipAC (2)

wobei M(d) als die Menge der Atomtupel (i,j) mit einer topologischen Distanz d zueinander

definiert ist.

Gleichung (2) lässt nun die Berechnung von Korrelationswerten für jede gegebene Distanz

und jede bekannte Atomeigenschaft zu. Es wird lediglich die Strukturformel des zu

kodierenden Moleküls benötigt. Betrachtet man zum Beispiel Distanzen von 0 bis 10

Bindungen und die Autokorrelation einer Atomeigenschaft, so hat der als Deskriptor

verwendete Korrelationsvektor folgende Datenstruktur:

AV = (AC0, AC1, AC2, AC3, AC4, AC5, AC6, AC7, AC8, AC9, AC10)

Der Deskriptor entspricht der Aneinanderreihung der Korrelationswerte zu den auftretenden

Atomeigenschaft bei einer Distanz d. Erfahrungsgemäß kann (2) aufgrund der

Aufsummierung zu besonders großen Korrelationswerten führen, wenn eine große Anzahl

von Atompaaren mit einer Distanz d vorliegt. Dies kann sich gerade bei dem Vergleich von

unterschiedlich großen Molekülen negativ auswirken, weshalb dies meist durch eine

Mittelung über die Anzahl der Atompaare eliminiert wird. Die erweiterte Formel lautet nun:

∑∈

⋅=)(,

)()(1

dMjid jpip

nAC , (3)

wobei n der Anzahl der Atompaare mit einer Distanz d entspricht.


© modlab® 2007

Beispiel

In der Tabelle 1 soll die Berechnung des Autokorrelationsvektors für zwei Moleküle A und B

beispielhaft gezeigt werden. Berechnet werden die Autokorrelationswerte für die

Elektronegativität der vorliegenden Atome (ENC=2,5; ENO = 3,5) mit einer Distanz von 1 bis

4 Bindungen unter Benutzung der Formel (3).

Tabelle 1: Berechnung der Korrelationswerte der Elektronegativität für zwei Moleküle A und B.

Verbindung A B

Molekülstruktur

OH

OH

OH

O

Name Butan-1,3-diol Buttersäure

Vereinfachte Struktur mit

willkürlichem Atomindex

O

O

1 23

45

6

O

O

12

34

5

6

AC1 7,25 7,25

AC2 7,75 8,45

AC3 4,75 4,75

AC4 3,7 2,5

Als Ergebnis erhält man folgende Vektoren:

VerbindungA 7,25 7,75 4,75 3,7

VerbindungB 7,25 8,45 4,75 2,5

Die Deskriptoren unterscheiden sich nur in den Funktionswerten für die Distanzen 2 und 4,

weil immer bei derjenigen Verbindung ein größerer Wert auftritt, bei dem die Sauerstoffe in

dem entsprechenden Abstand zueinander vorlagen. Eine Interpretation der resultierenden

Autokorrelationsvektoren ist in folgender Weise möglich: In der Verbindung A liegen mehr


© modlab® 2007

Atome mit der betrachteten Eigenschaft in einem Abstand von 4 Bindungen vor, als in der

Verbindung B. Hier tritt dies bei einer Distanz von 2 Bindungen auf.

Die betrachteten Eigenschaften müssen nicht direkt miteinander korreliert werden. In

bestimmten Fällen ist es ratsam modifizierte Werte der numerischen Eigenschaften zu

multiplizieren. Hier können unter anderem der Betrag, der Logarithmus u.ä. verwendet

werden.

Oft ist zudem die Verteilung von mehreren Atomeigenschaften auf den Molekülen zur

Beschreibung der biologischen Aktivität wichtig. Hier bedient man sich der

Kreuzkorrelation, bei der verschiedene Eigenschaften miteinander korreliert werden. Die

verwendete Formel lässt sich in folgender Weise zusammenfassen:

∑∑∈∈

⋅=)(,,

)()(1

dMjiPqpd jqip

nAC , (4)

wobei die Funktionen p und q den Elementen aus einer Funktionsmenge P mit allen

betrachteten Eigenschaften entspricht.

Autokorrelationsvektoren werden nicht nur zur Kodierung von zweidimensionalen

Molekülstrukturen verwendet. Eine Erweiterung der Methode zur Benutzung von 3-D

Molekülinformationen ist durch ein simples Austauschen der topologischen Distanz mit der

euklidischen Distanz möglich. Die Korrelationsvektoren kodieren dann für eine Beschreibung

der Eigenschaftsverteilungen im dreidimensionalen Raum, die durch das Molekül gegeben ist,

unter Berücksichtigung der Atomdistanzen in Ångstrom.

Anwendung von Korrelationsvektoren

Die Kodierung von Korrelationsvektor-basierten molekularen Deskriptoren liefert gleichgroße

Vektoren zur Beschreibung der vorliegenden Substanzen. Bei der Betrachtung von k

Atomeigenschaften und l Distanzen erhält man einen )( lk ⋅ -dimensionalen Deskriptor für

jedes Molekül. Zur Berechnung der Ähnlichkeit von zwei Molekülen bedient man sich

einfachen mathematischen Vektorfunktionen, wie der euklidischen Norm zweier Vektoren:

∑=

−=−n

iii baba

1

2)(vv

, (5)

wobei n der Dimensionalität des Deskriptors entspricht.


© modlab® 2007

Falls man mehrere Eigenschaften bei der Korrelation in Betracht zieht, können Fälle

auftreten, die sich aufgrund der Werte je Eigenschaft in (5) stärker auf die

Ähnlichkeitsanalysen auswirken. Ein einfaches Beispiel bildet zum Beispiel die Korrelation

der Elektronegativität von Atomen im Vergleich zu Korrelationen der Atommasse, welche

wesentlich höher ausfallen kann. Hier bietet es sich an, nach der Berechnung des

Korrelationsvektors auf eine Normierung des Vektors zurückzugreifen (z.B. eine Skalierung

der Werte von 0 bis 1). Neben Ähnlichkeitsanalysen werden Korrelationsvektor-basierte

Deskriptoren für viele weitere chemieinformatische Methoden verwendet.

Literatur

Moreau, G., Broto, P., 1980, Autocorrelation of Molecular Structures, Application to SAR Studies. Nouv. J. Chim., 4, 757. Wagener, M., Sadowski, J., Gasteiger, J., 1995, Autocorrelation of Molecular Surface Properties for Modelung Corticosteroid Binding Globulin and Cytosolic Ah Receptor Activity by Neural Networks. J. Am. Chem. Soc., 117, 7769. Schneider, G., Clément-Chomienne, O., Hilfiger, L., Schneider, P., Kirsch, S., Böhm, H.-J., Neidhart, W., 2000, Evolutionäres De-novo-Design bioaktiver Moleküle: ein Ansatz zum virtuellen Screening. Angew. Chem., 112, 4305. Schneider, G., Wrede, P., 1998, Artificial neural networks for computer-based molecular design. Prog. Biophys. Mol. Biol., 70, 175.

NIPALS-Algorithmus 29

© modlab® 2007

NIPALS-Algorithmus

Die Hauptkomponentenanalyse (Principal Component Analysis, im Folgenden mit „PCA“

abgekürzt) ist eine statistische Methode zur Analyse von multidimensionalen Daten. Mit Hilfe

dieser Technik ist es unter anderem möglich, die Dimensionalität eines Datensatzes mit einem

definiertem Informationsverlust zu reduzieren. Die Basis der PCA ist die Berechung von

Eigenvektoren. Diese haben als Eigenschaft, dass sie alle orthogonal zu einander sind. Sie

treten immer paarweise mit dem dazugehörigen Eigenwert auf, anhand welcher man die

Eigenvektoren sortieren kann. Der Eigenvektor mit dem größten Eigenwert ist als die erste

Hauptkomponente definiert. Mit ihr ist es möglich, die stärkste Beziehung in der

Datenmatrix aufzuzeigen. Aufgrund dieser Eigenschaften ist es möglich, die Daten einer

Matrix durch Ihre Projektion auf die Hauptkomponenten (Eigenvektoren mit den größten

Eigenwerten) darzustellen. Bei einer geschickten Auswahl von Eigenvektoren können so

dieselben Daten mit einer geringeren Dimensionalität wiedergegeben werden. In der Regel

wählt man hierzu alle Eigenvektoren mit einem Eigenwert > 1 aus. Zur Berechnung der

Eigenvektoren greifen viele Anwender auf die integrierten Lösungen der verwendeten

Programmpakete, wie Java oder Matlab, zu. Hier werden meist numerisch exakte,

mathematische Verfahren wie der QZ-Algorithmus, die Potenzmethode, die Inverse Iteration,

das Lanczos-, Arnoldi- oder Jacobi-Davidson-Verfahren verwendet. Alle numerischen

Verfahren haben gemeinsam, dass die Komplexität bei der Anwendung auf großen Matrizen

enorm steigt. Dadurch wird die Berechnung der Eigenvektoren sehr aufwändig und die

Methoden nicht praktikabel in der Anwendung.

Eine Alternative bilden approximative Lösungen. Ein populäres Beispiel zur

Eigenvektorabschätzung ist der 1966 von Herman Wold entwickelte NIPALS-Algorithmus

(nonlinear iterative partial least squares), welcher im Folgenden erklärt wird. Der NIPALS-


© modlab® 2007

Algorithmus wird nicht nur in der Hauptkomponentenanalyse angewendet.

Weiterentwicklungen und Derivate des Algorithmus finden sich auch in der multiplen

Regression, der Faktorzerlegung und der kanonischen Korrelationsanalyse wieder.

Prinzip und Schritte des Algorithmus

Die Hauptkomponenten werden nacheinander durch die iterative Anwendung der Schritte 1-9

(Tabelle 1) gewonnen.

Tabelle 1: Die Arbeitsschritte des NIPALS-Algorithmus.

Schritt Funktion Erklärung

1 ixs =: Wähle eine beliebige Spalte der Matrix X

als den initialen Scoringvektor s.

2

ss

sXl

T

T

= Projiziere X auf s, um den entsprechenden

Loadingvektor zu erhalten.

3

l

ll =

Normalisiere l auf die Länge 1.

4 ssold =: Zwischenspeichern des aktuellen Scores.

5

ll

lXs

T=

Projiziere X auf l, um den entsprechenden

neuen Scoringvektor zu erhalten.

6 ssd old −=: Bestimme die Differenz d des neuen und

alten Scores.

7 ?td < Falls d größer als ein vorher definierter

Schwellwert t (z.B. 10-6) ist, gehe zurück zu

Schritt 2. Ansonsten fahre fort bei Schritt 8.

8 TslXE −=: Entferne den Anteil der gefundenen

Hauptkomponente aus der Datenmatrix X.

9 EX =: Um weitere Hauptkomponenten zu finden,

wiederhole den Algorithmus ab Schritt 1 mit

der neuen Datenmatrix.


© modlab® 2007

Der NIPALS-Algorithmus approximiert eine Datenmatrix X durch die kleinere Scoringmatrix

S und die Loadingmatrix LT mit Hilfe der Matrixmultiplikation TSLX = (Abb. 1). Das

hochgestellte „T“ an einer Matrix verdeutlicht, dass die jeweilige Matrix transponiert wird.

Abbildung 2: Prinzip des NIPALS-Algorithmus. Ziel ist die Reproduktion der Datenmenge X durch eine Matrixmultiplikation von S und LT, die der Scoringmatrix und der transponierten Loadingmatrix entsprechen. Die Rechtecke skizzieren die jeweiligen Matrizen mit den an den Ecken angegebenen Dimensionen.

Literatur

Wold, H., Estimation of principal components and related models by iterative least squares, Multivariate Analysis (Ed., Krishnaiah, P. R.), Academic Press, NY, pp. 391-420 (1966). Wold S., 1974, A theoretical foundation of extrathermodynamic relationships (linear free energy rela-tionships), Chemica Scripta, 5, 97-106. Trygg, J., 2003, Chemometrics, Umeå University, Sweden, URL: http://www.chemometrics.se/editorial/sep2003.html#4

Graphentheorie 33

© modlab® 2007

Grundlagen Graphentheorie

Graphen

Der Graph als mathematisches Modell und Datenstruktur hat in der Cheminformatik eine

zentrale Bedeutung zur Repräsentation von chemischen. Seien FE die Menge der

Kanteneigenschaften und FV die Menge der Knoteneigenschaften.

Ein Graph ist ein 4-Tupel G=(V, E, µV, µE), wobei

• V eine endliche Menge der Knoten (engl. vertices),

• E ⊆ V×V die Menge der Kanten (engl. edges),

• µV: V→ FV eine Funktion ist, die Knoten Eigenschaften zuweist,

• µE: E→ FE eine Funktion ist, die Kanten Eigenschaften zuweist.

Man unterscheidet zwischen gerichteten und ungerichteten Graphen. Für ungerichtete

Graphen gilt die Bedingung:

• ∀ v,w ∈V: (v,w) ∈ E ⇔ (w,v) ∈ E

Die Kanten eines gerichteten Graphen haben also einen Anfangs- und einen Endpunkt, daher

werden sie häufig als Pfeile dargestellt (Abbildung ). Im ungerichteten Graphen ist es gleich,

welcher Knoten der Anfangs- und welcher der Endpunkt einer Kante ist.

Graphentheorie 34

© modlab® 2007

Abbildung 1: Beispiel a) für einen ungerichteten und b) für einen gerichteten Graphen.

Stellt ein Graph das Modell eines Moleküls dar, spricht man vom molekularen Graphen.

Dabei repräsentieren die Knoten die Atome und die Kanten die Bindungen zwischen den

Atomen. Die Knoten- und Kanteneigenschaften sind in diesem Fall Atom- und

Bindungseigenschaften, z.B. das Elementsymbol oder die Bindungsordnung. Die graphische

Repräsentation eines Moleküls ist in der organischen Chemie ein Standard (Abbildung ).

O

OH

Abbildung 2: Molekularer Graph von Tetrahydrocannabinol.

Da es meistens nicht von Bedeutung ist, welches Atom der Ausgangspunkt und welches der

Endpunkt einer Bindung ist, handelt es sich beim molekularen Graphen um einen

ungerichteten Graphen.

Es gibt mehrere Möglichkeiten, einen Graphen im Rechner zu repräsentieren:

• die Adjazenzliste

• die Adjazenzmatrix,

• und die Inzidenzmatrix.

Die Datenstrukturen mit den darauf operierenden effizienten Algorithmen werden in Java von

der Graphenbibliothek jgrapht (http://jgrapht.sourceforge.net/, Version 0.8.3) bereitgestellt.

Graphentheorie 35

© modlab® 2007

Subgraphen, Pfade, Zyklen, Cliquen

Sei G=(V, E, µv, µE) ein Graph. Dann heißt G’ = (V’, E’, µv’, µE’) ein Subgraph von G, wenn

• V’⊆V

• E’⊆E

• µV’(v) = µV(v) wenn v ∈ V’

• µE’(e) = µE(E) wenn e ∈ E’.

Ein Pfad in einem Graphen ist definiert als eine Reihe von Kanten e1, e2,..., en für die gilt:

• ei und ei+1 haben einen gemeinsamen Knoten.

• wenn für alle i ∉ {1, n} gilt: ei hat einen gemeinsamen Knoten mit ei-1 und den

anderen mit ei+1.

Ein Zykel (Kreis) ist ein Pfad, für den gilt

• e1 und en haben einen gemeinsamen Knoten.

ein ungerichteter Graph

ein Pfad (fett gekennzeichnet)

ein Kreis (fett gekennzeichnet)

eine Clique (fett gekennzeichnet)

Abbildung 3: Übersicht über die Definitionen zur Graphentheorie.

Ein Graph heißt vollständig, wenn jeder Knoten mit jedem anderen Knoten verbunden ist.

Eine Clique ist definiert als ein vollständiger Subgraph. Die Größe der Clique ist definiert als

die Anzahl der enthaltenen Knoten. Eine maximale Clique ist echt nicht in einer größeren

Graphentheorie 36

© modlab® 2007

Clique enthalten. Eine größtmögliche Clique ist eine Clique der Größe s, wobei gilt: der

Graph enthält keine weiteren Cliquen größer s. Alle Begriffe sind mit Beispielen in der

Abbildung zusammengefasst.

Graphentheoretische Probleme und ihre Lösungen in der

Cheminformatik

Subgraph-Isomorphismus

Die Suche von Substrukturen in Molekülen kann auf graphentheoretische Methoden

zurückgeführt werden. In diesem Fall ist sie komplett innerhalb eines anderen gut

untersuchten Problems enthalten: dem Problem des Subgraph-Isomorphismus.

Eine bijektive Funktion f: V→ V’ ist ein Graph-Isomorphismus vom Graphen G=(V, E, µV,

µE ) zum Graphen G’ = (V’, E’, µv’, µE’), wenn gilt:

• µV(v) = µV’(f(v)) für alle v∈V.

• für jede Kante e = (v,w) ∈ E existiert eine Kante e’ = (f(v),f(w)) ∈ E’, für die gilt µE(e)

=µE(e’), und für jede Kante e’ = (v’,w’) ∈ E’ existiert eine Kante e = (f-1(v’),f-1(w’)) ∈

E für die gilt µE(e’) =µE(e).

Eine injektive Funktion f: V→ V’ ist ein Subgraph-Isomorphismus von G zu G’ wenn gilt:

• es existiert ein Subgraph S ⊆ G’ so dass f ein Graph-Isomorphismus von G zu S ist.

Informell kann man das Problem des Subgraph-Isomorphismus so formulieren:

• „Gegeben seien die Graphen G1 und G2. Enthält G1 eine Kopie von G2 als Subgraph?

Können wir also eine Teilmenge der Knoten von G1 finden, die zusammen mit ihren

Kanten in G1 eine exakte Kopie von G2 bilden, wenn wir die Zuordnung zwischen den

Knoten und Kanten von G2 und den Knoten und Kanten des Teilgraphen von G1

entsprechend wählen?“

Graphentheorie 37

© modlab® 2007

Das Subgraph-Isomorphismus-Problem ist NP-vollständig. Folglich konnte kein Algorithmus

geschrieben werden, der den Subgraph-Isomorphismus in polynomialer Zeit findet. Dennoch

wurden Algorithmen gefunden, die eine akzeptable mittlere Laufzeit besitzen. In der

Chemieinformatik hat der Ullman-Algorithmus breiten Zuspruch gefunden.

Der Ullman-Algorithmus, der in 1976 von J.D. Ullman publiziert wurde, ist einer der

schnellsten Algorithmen für das Problem des Subgraph-Isomorphismus. Insbesondere eignet

er sich zum Auffinden von Subgraph-Isomorphismen auf Graphen, die in ihren topologischen

Eigenschaften den Molekulargraphen ähneln. Da der Original-Algorithmus für nicht

markierte Graphen geschrieben wurde, wird im Folgenden die Version für markierte Graphen

aus vorgestellt.

Der Ullman-Algorithmus vereinigt die Backtracking-Suche mit der „forward-checking“-

Technik. Der Pseudocode gibt Aufschluss über die genaue Funktionsweise des Ullman-

Algorithmus.

function ullman(G=(V,E,µV,µE), G’=(V’,E’,µ’V,µ’E))

//finds all subgraph isomorphisms of an arbitrary graph G in a

graph G’//

P=(pij): a |V|×|V’| matrix

F: set of vertex pairs, partial matching

01:for i = 1 to |V|

02: for j = 1 to |V’|

03: if µV(vi)= µ’V(v’j)

04: then pij = 1

05: else

06: pij = 0

07: end;

08:end;

09:backtrack(P,1,∅)

end function

Graphentheorie 38

© modlab® 2007

function backtrack(P,i,F)

10: if i > |V|

11: then F is an isomorphism

12: else

13: for all pij = 1 and j = 1 to |V’|

14: F ← F ∪ {vi,v’j};

15: P’ ← P;

16: for all k > i

17: p’ij = 0;

18: end;

19: if forwardchecking(P’,i,F) = true

20: then backtrack(P’,i+1,F);

21: F ← F\{vi,v’j};

22: end;

end function

function forwardchecking(P,i,F)

23: for k = i+1 to |V|

24: for l = 1 to |V’|

25: if pkl = 1

26: if consistent({vk,v’l},F)= false

27: pkl = 0

28: end;

29: end;

30: if there exists a row k in P such that pkl = 0 for l = 1..

|V’|

31: then return false

32: else

33: return true

end function

Graphentheorie 39

© modlab® 2007

function consistent({vk,v’l},F)

34: if for each {v,v’} ∈ F

35: if there is an edge (vk,v) ∈ E then there must be an

edge (v’,v’l) ∈ E’ with µE(vk,v)= µ’E(v’,v’l)

36: if there is an edge (v’,v’l) ∈ E’ then there must be

an edge (vk,v) ∈ E with µE(vk,v)= µ’E(v’,v’l)

37: then return true

38: else

39: return false

end function

Seien G=(V, E, µV, µE) und G’ = (V’, E’, µv’, µE’) zwei Graphen. Zunächst wird eine |V|×|V’|

Matrix P = (pij) erschaffen (Zeilen 1-8), für die gilt:

• wenn µV(vi)= µV’(v’i), dann ist pij= 1,

• sonst ist pij= 0.

In P steht also immer dann eine 1 in der i-ten Spalte und j-ten Zeile, wenn der Knoten vi den

knoten v’i matcht (Abbildung ).

Abbildung 4: a) Querygraph, b) Molekülgraph, c) daraus resultierende Matrix P.

Graphentheorie 40

© modlab® 2007

Danach wird die rekursive Prozedur backtrack aufgerufen (Zeile 09). Als Initialparameter

werden die Matrix P, Spaltenzähler i=1 und die Menge der erfolgreichen Matchings F=∅

übergeben.

Die Prozedur backtrack funktioniert wie folgt:

• ist der Spaltenzähler i > |V|, so ist F ein Isomorphismus (Zeilen 10 und 11). Dies ist

dann der Fall, wenn durch die rekursiven Aufrufe jedem Knoten aus V ein Knoten aus

V’ zugewiesen wurde.

• ist i ≤ |V|, dann ist der komplette Subgraph-Isomorphismus noch nicht gefunden. In

diesem Fall wird für jede 1 in der i-ten Spalte rekursiv backtrack aufgerufen. Als

Parameter werden der inkrementierte Spaltenzähler (i+1; Zeile 20), der bereits

gefundene Teil des Subgraphisomorphismus in Form der Menge F und die durch die

Prozedur forwardchecking bearbeitete Matrix P’.

Die Prozedur forwardchecking minimiert die Anzahl der Rekursionsschritte. Dazu wird

jede 1 in der Matrix P, also jedes „matchende“ Knotenpaar (v,v’) auf Konformität mit jedem

Element aus der Menge der bereits gefundenen „matchenden“ Knotenpaare F überprüft. Zwei

Knotenpaare (v,v’) und (w,w’) sind konform, wenn gilt:

• wenn es eine Kante (v,w) gibt, muss es auch eine Kante (v’,w’) geben für die gilt:

µE(v,w) = µ’E(v’,w’), und

• wenn es eine Kante (v’,w’) gibt, muss es auch eine Kante (v,w) geben für die gilt:

µE(v,w) = µ’E(v’,w’).

Somit erfolgt in der Prozedur forwardchecking nicht nur die Minimierung der

Backtrackschritte, sondern auch die Überprüfung, ob F auch tatsächlich ein Subgraph-

Isomorphismus von G zu G’ ist. Außerdem wird die Rekursion unterbrochen, wenn in P eine

Zeile vorhanden ist, die nur 0er enthält, also kein Subgraph-Isomorphismus in diesem Ast des

Rekursionsbaums gefunden werden kann.

Maximum Common Substructure

Eine weitere graphentheoretische Anwendung in der Cheminformatik ist die Berechnung der

Isomorphie größtmöglicher gemeinsamer Teilgraphen (Maximum Common Substructure,

MCS). Die MCS zweier Molekülgraphen ist in Abbildung 5 dargestellt.

Graphentheorie 41

© modlab® 2007

O

HO

N

OHN

N

OH

Abbildung 5: MCS zweier Molekülgraphen

Die MCS kann zur Berechnung der Ähnlichkeit zweier Moleküle benutzt werden und ist wie

folgt definiert:

Seien G1 und G2 zwei Graphen, dann ist MCS eine Struktur (S1, S2), so dass gilt

• S1 ist Subgraph von G1

• S2 ist Subgraph von G2

• S1 und S2 sind isomorph

• es gibt kein S1’ für das gilt: S1 ist Subgraph von S1’

• es gibt kein S2’ für das gilt: S2 ist Subgraph von S2’

Um alle MCS zweier Graphen aufzuzählen, bedient man sich u.a. der Cliquendetektion auf

dem Assoziationsgraphen.

Cliquendetektion auf dem Assoziationsgraphen

Die Lösung des MCS-Problems durch Cliquendetektion auf dem Assoziationsgraphen besteht

im wesentlichen aus zwei Schritten:

• Erstellung des Assoziationsgraphen,

• Cliquendetektion auf dem Assoziationsgraphen.

Für zwei Graphen G=(V, E, µV, µE) und G’ = (V’, E’, µv’, µE’) ist der Assoziationsgraph

(Abbildung 6) definiert als ein unmarkierter Graph GA = (VA, EA) mit VA ⊆ V × V’ und EA ⊆

VA × VA. vA ∈ VA heißt a-Knoten (Assoziationsknoten), eA ∈ EA heißt a-Kante

(Assoziationskante) von GA. Für GA muß gelten:

Graphentheorie 42

© modlab® 2007

• VA = {(v,v’)| v ∈ V, v’ ∈ V’, µV(v) = µV’(v’) }

• EA besteht aus a-Kanten eA= (vA,wA) mit vA=(v,v’), wA=(w,w’), so dass gilt

o v ≠ w und v’ ≠ w’

o wenn eine Kante e = (v,w) ∈ E existiert, muss eine Kante e’ = (v’,w’) ∈ E’

existieren und µE (e) = µE’(e’) muss gelten.

o wenn keine Kante e = (v,w) ∈ E existiert, darf auch keine Kante e’ = (v’,w’) ∈

E’ existieren.

Abbildung 6: a) Querygraph. b) Molekülgraph. c) Assoziationsgraph aus a) und b), grau eingezeichnet: Clique der Größe 4.

Wenn also die a-Knoten (v1,v’1), (v2,v’2)… (vn,v’n) für 1≤ n ≤ |V|⋅|V’| im Assoziationsgraphen

paarweise miteinander verbunden sind, ist der Subgraph in V, der durch die Knoten v1, v2,…,vn

induziert ist, isomorph zum Subgraphen in V’, der durch die Knoten v’1, v’2,…,v’n induziert

ist. Der Isomorphismus ist durch (v1,v’1), (v2,v’2)… (vn,v’n) explizit gegeben.

Abbildung 6 zeigt die Erstellung eines Assoziationsgraphen aus dem Molekülgraphen g1

(Alanin) und einem Querygraphen g2 (Carbonylfunktion verbunden mit einem Stickstoff über

eine Kohlenstoffbrücke). Die Atomtypen entsprechen dabei den Knoteneigenschaften, die

Bindungsordnungen den Kanteneigenschaften. Nun wird z.B. ein a-Knoten (3,1’) erstellt, weil

der Knoten 3 dieselbe Eigenschaft trägt wie der Knoten 1’, nämlich O (Sauerstoff). Einen a-

Knoten (3,4’) gibt es nicht, da der Knoten 3 die Eigenschaft O (Sauerstoff) und der Knoten 4’

Graphentheorie 43

© modlab® 2007

die Eigenschaft N (Stickstoff) trägt. Ähnlich verfährt man bei der Erstellung der a-Kanten.

Die a-Knoten (3,1’) und (2,2’) werden durch eine a-Kante verbunden, da sowohl in g1 die

Knoten 3 und 2 als auch in g2 die Knoten 1’ und 2’ durch eine Kante mit der Eigenschaft

„Doppelbindung“ verbunden sind. Die a-Knoten (3,1’) und (6,4’) sind verbunden, weil es in

g1 keine Kante (3,6) und in g2 keine Kante (1’,4’) existiert. Der komplette a-Graph enthält

eine Clique der Größe 4, nämlich {(3,1’), (4,3’),(2,2’),(6,4’)} (in Abbildung , rot

hervorgehoben). Offensichtlich entspricht sie dem MCS von g2 zu g1.

Der zweite Schritt besteht in der Suche der Cliquen auf dem a-Graphen. In der Literatur gibt

es vielfältige Ansätze zur Lösung dieses Problems. Der schnellste vollkombinatorische Ansatz

ist der Bron-Kerbosch-Algorithmus. Es handelt sich dabei um einen rekursiven Backtracking-

Algorithmus, der alle maximalen Cliquen in einem Graphen aufzählt. Im Wesentlichen

operiert der BK-Algorithmus auf drei Mengen C, P und S.

Die Menge C (Clique) ist die Menge der Knoten, die bereits als komplett verbunden gefunden

wurden, also eine Clique bilden. Die Menge P (Potential) enthält diejenigen Knoten, die zu

dem zuletzt zu C hinzugekommenen Knoten adjazent sind, folglich also für die Erweiterung

der Clique benutzt werden können. Die Menge S enthält all die Knoten, die nicht mehr für die

Erweiterung von C benutzt werden können, weil alle Cliquen, die diese Knoten enthalten,

bereits gefunden wurden. Der Pseudocode gibt die Funktionsweise des BK-Algorithmus

wieder.

Der Algorithmus wird mit C,S = ∅ und P=V initialisiert. Wenn die Mengen S und P leer sind,

ist C eine maximale Clique und kann ausgegeben oder gespeichert werden (Zeile 03). Ist dies

nicht der Fall, so wird jeder Knoten aus P in der Schleife (Zeilen 04-11) bearbeitet. Ein

ausgesuchter Knoten ui wird aus P entfernt (Zeile 05), P wird nach P’ (Zeile 06) und S nach

S’ (Zeile 07) kopiert, um für den rekursiven Aufruf verwendet zu werden. Knoten, die mit ui

verbundenen sind, werden in N gespeichert (Zeile 08). Der rekursive Aufruf erfolgt mit der

Menge C∪{ui} (die bereits gefundene Clique wird um ui erweitert), P∩N (nur die Knoten in

P, die mit ui verbunden sind, dürfen für die potentielle Erweiterung der Clique verwendet

werden) und S∩N (alle Knoten, die bereits betrachtet wurden und mit ui verbunden sind)

(Zeile 09). Anschließend wird ui der Menge S zugefügt (Zeile 10).

Graphentheorie 44

© modlab® 2007

function find_all_cliques(C,P,S)

//enumerates all cliques in an arbitrary graph G //

C: set of vertices belonging to the current clique

P: set of vertices which can be added to C

S: set of vertices which are not allowed to be added to C

N[u]: set of vertices which are adjacent to vertex u in G

01: Let P be the set {u1,...,uk}⊆ V;

02: if P = ∅ and S = ∅

03: then C is a maximal clique

04: else for i ← 1 to k

05: do P ← P\{ui};

06: P’ ← P;

07: S’ ← S;

08: N ← {v∈V | {ui,v}∈E};

09: find_all_cliques(C∪{ui},P’∩N,S’∩N);

10: S ← S ∪ {ui};

11: end;

12:end;

end function

Der BK-Algorithmus wurde bereits in vielen bioinformatischen Zusammenhängen verwendet

und modifiziert. Die besondere Geschwindigkeit der sog. „Version2“ des BK-Algorithmus

liegt in einer Heuristik bei der Wahl des Knotens ui aus P. Diese motiviert sich aus der

Verwendung der Menge S, da man die Abbruchbedingung für die Ausgabe der Clique auch so

formulieren kann: Es gibt einen Knoten s in S, der mit allen Knoten in P verbunden ist. Ein

solcher Knoten s kann niemals aus S entfernt werden. Somit haben wir eine Möglichkeit,

vorzeitig vorherzusagen, ob ein Ast im Backtrackingbaum zum Erfolg führt oder nicht.

Allgemein nennt man diese Backtracking-Technik die „branch-and-bound“ Methode, das

Abbruchkriterium ist das „branch-and-bound“ Kriterium. Die „Version2“ des BK-

Algorithmus sucht ui aus P so aus, dass das „branch-and-bound“ Kriterium möglichst schnell

eintritt. Dies geschieht in folgender Weise: alle Knoten in S sind mit einem Zähler assoziiert,

welcher die Anzahl der Kanten zu den Knoten in P wiedergibt. Verschiebt man nun einen

Knoten p aus P nach S, so vergrößern sich die Zähler aller Knoten in S höchstens um 1 und p

bekommt einen neuen Zähler. Das „branch-and-bound“ Kriterium kann nun folgendermaßen

Graphentheorie 45

© modlab® 2007

formuliert werden: Es gibt einen Knoten s in S mit Zähler z =|P|. Bron und Kerbosch zeigen,

dass man dieses Kriterium am schnellsten erreicht, wenn man ui aus P so auswählt, dass ui mit

demjenigen Knoten s in S durch eine Kante verbunden ist, welcher den höchsten Zähler z hat.

Dadurch erhöht man in jedem Schritt den Zähler von s um 1, das „branch-and-bound“

Kriterium wird auf dem schnellst möglichen Wege erreicht.

Literatur

Graphentheorie, Komplexizität, grundlegende Algorithmen G. Valiente (2002), Algorithms on Trees and Graphs, Springer Berlin Heidelberg. Ottmann T, Widmayer P (1996). Algorithmen und Datenstrukturen, Spektrum Akademischer Verlag GmbH, Heidelberg, 537ff. Garey R, Johnson DS (1979). Comuters and Intractability: a Guide to the Theory of NP-Completeness, W. H. Freeman and Company, San Francisco. Schöning U (2001), Algorithmik, Spektrum Akademischer Verlag GmbH, Heidelberg. Algorithmen zur Substruktursuche Barnard J (1993). "Substructure Searching Methods - Old and New." J. Chem. Inf. Comp. Sci. 33, 532-538. Messmer B (1995). Efficient Graph Matching Algorithms, Universität Bern, Bern, 11ff. Ullmann J (1976). "Algorithm for Subgraph Isomorphism." Journal of the ACM 23, 31-42. Cliquendetectionsalgorithmen und MCS Bron C, Kerbosch J (1973). "Agorithm 457: Finding all Cliques of an Undirected Graph." Communications of the ACM 16, 575-577. Koch I (2001). "Enumerating all connected maximal common subgraphs in two graphs." Theoretical Computer Science 250, 1-30. Raymond J, Gardiner E, Willett P (2002). "RASCAL: Calculation of graph similarity using maximum common edge subgraphs." Computer Journal 45, 631-644.

Graphentheorie 46

© modlab® 2007

Künstliche neuronale Netze 47

© modlab® 2007

Künstliche neuronale Netze

Künstliche neuronale Netze (engl.: artificial neural net, ANN) bilden eine wichtige Kategorie

im Bereich des maschinellen Lernens und der künstlichen Intelligenz. Ihre Aufgabe ist es,

anhand von Trainingsbeispielen eine Funktion (allgemein: ein Muster) zu erlernen und

anschließend das Erlernte auf unbekannte Beispiele anzuwenden. Ihre generelle Fähigkeit zur

Mustererkennung und Klassifizierung macht ANNs in großer Breite in der Praxis einsetzbar:

Von der Spracherkennung über Komprimierung und Rauschunterdrückung in der

Bildverarbeitung bis hin zur Vorhersage von Aktienkursverläufen wird auf künstliche

neuronale Netze zurückgegriffen.

ANNs sind in ihrem Aufbau und in ihrer Funktion biologischen neuronalen Netzen

nachempfunden. Das Zusammenspiel relativ simpler Untereinheiten, die in Anlehnung an die

biologische Terminologie als Neuronen bezeichnet werden, ergibt durch deren Vernetzung ein

komplexes System zur Informationsverarbeitung.

Die Art der Vernetzung der Neuronen ist ein Merkmal zur Unterteilung in verschiedene

Klassen künstlicher neuronaler Netze. So unterscheidet man im Wesentlichen zwischen

vollständig und teilweise verbundenen Netzen. Als Vertreter der vollständig verbundenen

Netze sei hier das Hopfield-Neural-Network genannt. In ihm ist jedes Neuron mit allen

anderen Neuronen verbunden.

Dieser Teil des Skriptes konzentriert sich hauptsächlich auf eine bestimmte Klasse der

teilweise verbundenen Netzwerke – den Feed-Forward-Netzen.

Genereller Aufbau von Feed-Forward-Netzen

Feed-Forward-Netze folgen einem schichtartigen Aufbau. Die oberste Schicht wird als Input-

Layer bezeichnet. Darauf folgt eine (theoretisch) beliebige Anzahl von sog. Hidden-Layer.

Die abschließende Schicht wird als Output-Layer bezeichnet. Abbildung 1 zeigt den


© modlab® 2007

allgemeinen Aufbau eines vierlagigen Feed-Forward-Netzes. Jede Schicht wird also bei der

Nomenklatur gezählt.

Abbildung 1: Aufbau eines vierlagigen Feed-Forward-Netzes.

Die namensgebende Eigenschaft der Feed-Forward-Netze ist, dass die

Informationsverarbeitung nur in eine Richtung (vom Input-Layer in Richtung des Output-

Layer) und zwischen direkt benachbarten Schichten erfolgt. Daraus resultiert die in

Abbildung 1 gezeigte Vernetzung zwischen den benachbarten Schichten: Jedes Neuron eines

Hidden-Layers ist mit allen Neuronen der darüber liegenden Schicht und mit allen Neuronen

der darunter liegenden Schicht verbunden.

Im Folgenden soll zunächst auf die einfachste Feed-Forward-Architektur - das Perzeptron -

eingegangen werden, um anschließend komplexere Beispiele zu betrachten.

Das Perzeptron

Das Single Layer Perzeptron

Das Single Layer Perzeptron, kurz SLP, stellt die einfachste Art eines Feed-Forward-Netzes

dar. Entwickelt wurde es 1958 von Frank Rosenblatt. Es besteht aus einem Input-Layer,

dessen Neuronen direkt mit einem Output-Neuron verbunden sind (Abbildung 2).

Zur Vermeidung von Missverständnissen: Nach obiger Definition handelt es sich beim SLP

um ein zweilagiges Netz, welches jedoch nur eine aktive Schicht besitzt und deshalb als

Single Layer bezeichnet wird.


© modlab® 2007

Die Verbindungen zwischen den Neuronen sind mit Gewichten versehen. Diese beschreiben

die Intensität, mit der der jeweilige Input auf das Ergebnis des Outputs einwirkt. In dem

Ouput-Neuron werden die Eingaben verarbeitet. Es wird eine gewichtete Summe, nämlich die

Summe der Produkte der Gewichte und der Inputwerte, berechnet. Weiterhin wird ein

zusätzlicher Wert, das sog. Bias, zur Summe hinzuaddiert. Das Bias-Neuron ist ein

zusätzliches Neuron, welches stets den Wert 1 ausgibt. Wie alle Neuronen besitzt es ein

Gewicht. In Abbildung 2 ist dieses mit Ө bezeichnet.

Abbildung 2: Schematische Darstellung eines Perzeptrons

Das Ergebnis dieser Berechnung wird in der sog. Aktivierungs- oder

Schwellenwertfunktion verarbeitet. Der Wert, den diese Funktion annimmt, bezeichnet den

Zustand des Neurons. Für eine binäre Schwellenwertfunktion gibt er an, ob ein Neuron

aktiviert oder deaktiviert ist (Abbildung 3), bzw. bei anderen Aktivierungsfunktionen, wie

stark die Aktivierung ist (Abbildung 6). Im Falle des SLPs, das mit einer binären

Schwellenwertfunktion arbeitet, gibt das Neuron als Output den Wert 1 aus, wenn das Neuron

aktiviert ist. Sollte das Ergebnis der gewichteten Summe unter dem gegebenen Schwellenwert

liegen, so geht das Neuron in einen deaktivierten Zustand über. Die Ausgabe beträgt in

diesem Fall -1. Der Schwellenwert liegt bei der in Abbildung 3 gezeigten binären

Schwellenwertfunktion bei 0.

θ+

= ∑

=

n

iii wx

1


© modlab® 2007

Abbildung 3: Binäre Schwellenwertfunktion des SLP. Der Schwellenwert der Funktion ist hier 0.

Mit Hilfe des Bias ist eine implizite Verschiebung des Schwellenwertes möglich. Wie wir

später sehen werden, ist das Training eines neuronalen Netzes durch die Veränderung der

Gewichte charakterisiert. Die Veränderung des Bias zieht eine Veränderung des

Schwellenwertes des Neurons nach sich. Sinkt das Bias-Gewicht, wirkt sich dies durch ein

Sinken der gewichteten Summe aus. Folglich erhöht sich dadurch implizit der Schwellenwert

für das Neuron, d.h. Eingaben und Gewichte müssen größer sein, um mit gleicher Stärke auf

die Aktivität des Neurons einzuwirken. Der Vorteil, die Veränderung des Schwellenwertes

über das Bias darzustellen, liegt darin, dass in Netzen mit mehr als einem Neuron für alle

Neuronen die gleiche Aktivierungsfunktion benutzt werden kann. Explizit bleibt der

Schwellenwert der Aktivierungsfunktion unverändert für alle Neuronen gleich. Durch das

Verändern des Bias-Gewichts während des Trainings kann er jedoch implizit für jedes Neuron

individuell verschoben werden. Ohne das zusätzliche Bias-Neuron benötigte man dafür für

jedes Neuron eine individuelle Aktivierungsfunktion, die zusätzlich während des Trainings

verändert werden müsste.

SLPs sind in der Lage, Klassifizierungsprobleme zu lösen, die linear separierbar sind. Die

Klassifizierung erfolgt hierbei implizit durch das Legen einer trennenden Hyperebene durch

die Punktmenge. Diese Hyperebene ist in ihrer Dimension stets um eins geringer als die

Dimension der zu erlernenden Funktion.

Explizit erfolgt die Trennung durch die vom Perzeptron erlernte Funktion, die orthogonal zu

dieser trennenden Hyperebene steht.


© modlab® 2007

Abbildung 4 illustriert die lineare Trennbarkeit. Hierbei ist die eingezeichnete Gerade die

trennende Hyperebene. Wie man sieht ist die zu erlernende Funktion zweidimensional,

während die trennende Hyperebene eine Gerade ist, die sich durch die allgemeine

eindimensionale Geradengleichung darstellen lässt. Sie ist in (1) gezeigt.

baxxf +=)( (1)

Ist das Klassifizierungsproblem nicht linear separierbar, kann es durch das SLP nicht gelöst

werden. Als populäres Beispiel kann hier die XOR-Funktion angeführt werden. Es kann durch

ein einzelnes Perzeptron nicht gelöst werden.

Abbildung 4: Ein linear separierbares Klassifizierungsproblem im zweidimensionalen Deskriptorraum.

Die Lösung des XOR-Problems

Für logische Operatoren wie AND oder OR ist das SLP in der Lage, durch das Legen einer

Geraden, die Punkte so zu separieren, dass eine Lösung möglich ist. Betrachtet man aber

Abbildung 5, so sieht man, dass dies für den XOR-Fall nicht möglich ist. Egal, wie eine

Gerade durch den Raum gelegt wird, ist es nicht möglich, dass diese die Punkte {(0,0),(1,1)}

von den Punkten {(0,1),(1,0)} trennt. Gezeigt wurde dies 1969 von Minsky und Papert.

Die Separierung des Raumes durch zwei Geraden würde zu der Lösung des Problems führen.

Dies ist allerdings nur mit einem mehrschichtigen Perzeptron lösbar.


© modlab® 2007

Abbildung 5: Das XOR-Problem. Die umrahmten Eingaben sollen in eine Klasse zusammengefasst werden, die anderen beiden in eine zweite Klasse.

Mehrlagige Netze

Ein einzelnes Neuron erstellt aus seinen Inputwerten und den zugehörigen Gewichten eine

gewichtete Summe, also eine lineare Funktion. Das SLP, welches im Wesentlichen ein

Neuron darstellt, kann also nur linear trennbare Klassifizierungsprobleme lösen. Wie am

Beispiel der XOR-Funktion gesehen, sind für das Perzeptron damit viele (sogar der

überwiegende Teil) der Klassifizierungsprobleme nicht lösbar. Schon das XOR-Problem

zeigt, dass der Übergang vom SLP zum Multi Layer Perzeptron (abgekürzt MLP) die Menge

der lösbaren Probleme vergrößert: MLPs sind generelle Funktionsapproximatoren.

Der Übergang zu mehrlagigen Netzen mit mehreren Neuronen pro Hidden-Schicht ergibt eine

Kombination linearer Gleichungen. Damit können beliebige Funktionen approximiert werden.

MLPs enthalten nichtlineare Aktivierungsfunktionen, z.B. die Sigmoidalfunktion und der

hyperbolische Tangens. Neben ihrem nichtlinearen Verlauf zeichnet diese Funktionen aus,

dass sie stetig und differenzierbar sind. Dies ist, wie wir noch sehen werden, für das Netz-

Training wichtig. Weiterhin bilden sie alle möglichen Ergebnisse der gewichteten Summe auf

das Intervall ]0,1[ bzw. ]-1,1[ ab. Abbildung 6 zeigt die entsprechenden Funktionsverläufe.

Achtung: Die Fähigkeit eines MLPs beliebige, auch nichtlineare Funktionen zu

approximieren, ist nicht durch die Wahl der Aktivierungsfunktion begründet, sondern durch

die Architektur des Netzes.


© modlab® 2007

Abbildung 6: Links die Sigmoidalfunktion, rechts der hyperbolische Tangens

Die Neuronen der Hidden-Layer geben also nicht direkt das Ergebnis der gewichteten Summe

an die nächste Schicht weiter, sondern ein zuvor durch die Aktivierungsfunktion

transformierten Wert.

Eine zentrale Frage beim Einsatz von ANNs zu Klassifizierung ist die der zu verwendenden

Netztopologie. Wie viele verdeckte Schichten werden benötigt? Wie viele Neuronen sind in

der jeweiligen Schicht optimal? Es gibt kein deterministisches mathematisches Verfahren, das

die optimale Netztopologie für ein gegebenes Problem bestimmt. Bei der Wahl der

Netztopologie stehen sich wie so oft zwei gegensätzliche Anforderungen gegenüber:

Einerseits soll das Netz die Möglichkeit haben, sich gut an die Trainingsbeispiele anzupassen.

Hierfür ist ein großes Netz von Vorteil. Andererseits soll das trainierte Netz möglichst gut

generalisieren, also ihm unbekannte Eingaben korrekt klassifizieren, und nicht einfach nur

den Trainingsdatensatz auswendig lernen. Diesem als Overfitting bezeichneten Problem

entgeht ein kompaktes Netz eher. Deshalb bleibt hier oft nur der trial-and-error-Ansatz.

Dennoch lässt sich generell sagen, dass in aller Regel mehr als zwei Hidden-Layer keinen

Sinn machen. Bereits Netze mit dieser Anzahl an verdeckten Schichten können potenziell alle

in der Praxis auftretenden Trennungsprobleme lösen. Eine weitere Vergrößerung der Netze

würde also keinen Vorteil für deren Mächtigkeit bringen, gleichzeitig aber die Gefahr des

Overfittings erhöhen und das Training verlangsamen.

Netztraining

Das Training eines ANN ist nichts anderes als die Anpassung der Gewichte zwischen den

Neuronen. Ziel ist es, die Gewichte so anzupassen, dass der Fehler des Netzes minimal wird.

Mit Fehler ist hier die Summe der Quadrate der Differenzen zwischen den vom Netz

berechneten Werten für die Trainingsmuster und deren gegebenen Zielwerten gemeint. Die


© modlab® 2007

Existenz der Zielwerte ist es, die das Training von ANNs zu einem Vertreter des überwachten

Lernens macht. Zusammengefasst führt dies zu Gleichung (2), wobei E den Fehlerwert des

Netzes summiert über alle k Trainingsmuster, ti den Zielwert des Trainingsmusters i und oi

dessen vom Netz berechneten Output-Wert darstellt.

∑=

−=k

iii towE

1

2)()(v

à min. (2)

Dies ist die während des Trainings zu minimierende Funktion. Macht man sich klar, dass in

die Berechnung der Output-Werte oi die Netzgewichte maßgeblich einfließen, ist leicht

ersichtlich, dass E eine Funktion der Gewichtswerte des Netzes ist (in (2) ist wv

der Vektor

aller Gewichtswerte des Netzes). Die Anzahl der Gewichte gibt also die Dimension des

Raumes an, in dem die Fehlerfunktion definiert ist. Abbildung 7 zeigt eine mögliche

Fehlerfunktion eines Netzes mit zwei Gewichten.

Auf die Gefahren in sehr hochdimensionalen Räumen und die unter anderem daraus folgende

Bevorzugung möglichst kompakter Netze soll später noch eingegangen werden.

Der populärste Trainingsalgorithmus für ANNs - der Backpropagtion of Errors (kurz

Backprop) - basiert auf der Idee des Gradientenabstieges. Dabei wird stets in Richtung des

steilsten Gefälles des Funktionsgraphen gelaufen. Bildlich kann man sich dies als einen

zufällig im Gebirge abgesetzten Skifahrer (das entspricht der zufällige Gewichtsinitialisierung

beim Erzeugen des Netzes) vorstellen, der anschließend eine Abfahrt immer entlang des

stärksten Gefälles startet, bis er im nächstgelegenen Tal landet.

Abbildung 7: Fehlerlandschaft, Quelle: [1]


© modlab® 2007

Der Bewegung des Skifahrers entspricht in unserem Problem eine Veränderung der Gewichte,

sodass sich der Fehler des Netzes verringert. Das Gradientenabstiegverfahren ermöglicht es

nun, für jedes Gewicht (jede Dimension) die Richtung (Vorzeichen) und die Intensität

(numerischer Wert) der Veränderung zu finden, um den Fehler zu verkleinern. Mathematisch

wird dies durch eine partielle Ableitung der Fehlerfunktion nach jedem Gewicht erreicht. Der

dabei entstehende Vektor von partiellen Ableitungen wird als Gradient oder Nabla-Operator

bezeichnet und ist in (3) definiert.

∂

∂

∂

∂

∂

∂=

nxxxg ,...,,

21

(3)

Backpropagation nimmt nun die Modifizierung der Gewichte proportional zum Gradienten

vor, wobei man den Proportionalitätfaktor als Lernrate η bezeichnet. Sie dient als Faktor der

linearen Modulierung des Einflusses des Gradienten. In unserem Bild des Skifahrers wäre

dies so etwas wie die Qualität seiner Skier: Je besser das Material (je höher η), desto schneller

fährt er bei gleichem Gefälle. Zusammengefasst ist dies in (4).

gw ⋅−=∆ ηv

(4)

wv

∆ stellt die Veränderung der einzelnen Gewichte zusammengefasst als Vektor dar, während

g den Gradienten der Fehlerfunktion E bezeichnet. Am hypothetischen Beispiel eines Netzes

mit nur einem Gewichtswert soll das Vorgehen des Gradientenabstieges verdeutlicht werden

(Abbildung 8).

Der Gradient ist die erste Ableitung der gezeigten Funktion nach dem Gewicht w. Hier wird

auch deutlich, warum in (4) der negative Gradient zur Gewichtsänderung herangezogen wird:

Der Funktionswert der ersten Ableitung ist bei w1 negativ, die Änderung des Gewichtes muss

jedoch positiv sein, um sich dem Minimum der Funktion zu nähern. Genau umgekehrt liegt

der Fall bei w2. Außerdem ist in Abbildung 8 ersichtlich, dass sich die Neigung der Geraden,

die ja die partielle Ableitung der Funktion an einer Position darstellt, auf die Stärke der

Gewichtsänderung auswirkt. Je stärker die Steigung (je höher der Wert der Ableitung), desto

größer ist auch die Veränderung des Gewichtes hin zum Optimum.


© modlab® 2007

Abbildung 8: Gradientenabstieg. w1 und w2 stellen zwei unterschiedliche Werte des Gewichtes w dar. Die Steigungen der Geraden entsprechen den Ableitungen der Funktion E(w) an diesen Positionen dar. Das zu findende Minimum der Funktion liegt bei w = 0.

Backpropagation stellt ein iteratives Verfahren zur Verfügung, um einen Gradientenabstieg

auf der Fehlerfunktion zu ermöglichen. Die Veränderung der Gewichtwerte wird dabei wie in

(5) bestimmt.

ijoldij

newij www ∆+= . (5)

�wij stellt die Veränderung des Gewichtwertes wij in einem Trainingsschritt dar. Die

Berechnung von �wij ist in (6) gezeigt und geht auf (4) und damit auf den Gradientenabstieg

zurück, wobei hier die Fehlerfunktion E partiell nach dem zu aktualisierenden Gewicht wij

abgeleitet wird. (Zur Wiederholung: Dies entspricht einer Komponente des Gradienten g aus

(4), der ja den Vektor aller partiellen Ableitungen darstellt).

ij

ij

ij xw

Ew ηδη =

∂

∂−=∆ . (6)

wij stellt hier das Gewicht zwischen Neuron i und Neuron j dar; xi ist der Ausgabewert des

Neurons i. Entscheidend für den Backpropagation-Algorithmus ist nun, dass sich diese

partiellen Ableitungen der Fehlerfunktion nach den Gewichten und damit δj effizient

berechnen lassen. Dabei muss allerdings unterschieden werden, ob es sich bei dem zu


© modlab® 2007

verändernden Gewicht um ein Gewicht zwischen dem untersten Hidden-Layer und der

Output-Layer handelt, oder um ein Gewicht zwischen höher liegenden Schichten. Die

Berechnung von δj unterscheidet sich also je nach dem, ob Neuron j ein Output-Neuron oder

ein Hidden-Neuron ist. Zusammengefasst ist dies in (7).

−

= ∑≤≤ nk

jkkj

jjj

j istonHiddenneurjfallswnetact

istonOutputneurjfallsotnetact

1

)('

))(('

δδ (7)

Hierbei entspricht tj im oberen Teil von (7) dem Zielwert und n im unteren Teil der Anzahl

der Neuronen in der Schicht unterhalb des Neurons j. Abbildung 9 soll der besseren Übersicht

und der Klärung der Variablenbezeichnungen dienen.

Abbildung 9: Kantengewichte sind rechts von ihrer Kante geschrieben, Ein- bzw. Ausgabewerte der Neuronen links davon. Die Bias-Neuronen sind durch „B“ gekennzeichnet, die Kästchen stellen die Aktivierungsfunktion dar.

Exemplarisch sollen die Berechnungen der �w-Werte für die Gewichte wij und vjk aus

Abbildung 9 demonstriert werden. Nehmen wir an, das angegebene Netz benutzt die in (8)

gezeigte Sigmoidalfunktion zur Aktivierung. Da wir nach (7) deren erste Ableitung benötigen,


© modlab® 2007

wird auch ersichtlich, warum eine Forderung an Aktivierungsfunktionen ihre

Differenzierbarkeit ist. Die erste Ableitung der Sigmoidalfunktion sig(x) ist, wie in (9)

ersichtlich, leicht zu errechnen.

xexsig

−+=

1

1)( . (8)

))(1()()(' xsigxsigxsig −⋅= . (9)

Damit ergibt sich nach (6) für ∆vjk (k ist Output-Neuron):

jjk o

k

knetsig

knetsig

ko

ktv ⋅−⋅⋅−⋅=∆

444444 3444444 21δ

η ))(1()()( ,

und für ∆ wij (j ist Hidden-Neuron):

iij inp

j

jnetsig

jnetsig

jkv

kw ⋅−⋅⋅⋅=∆

444444 3444444 21δ

δη ))(1()()( .

Aus (7) ist auch ersichtlich, dass zur Berechnung der δ-Werte der Neuronen einer Schicht die

δ-Werte der Neuronen der darunter liegenden Schicht benötigt werden. Einzig die δ-Werte

des Output-Layers sind direkt berechenbar. Der Netzfehler wird also durch den Algorithmus

vom Output-Layer über die Hidden-Layer von unten nach oben durchgereicht

(„backpropagation“). Dies gab dem Trainingsalgorithmus seinen Namen. Die Gewichte der

Bias-Neuronen werden wie alle anderen Gewichte nach dem beschriebenen Schema während

des Trainings angepasst.

Für Details der mathematischen Herleitung sei auf [1] verwiesen. Dort wird genau

beschrieben, wie die partielle Ableitung von E nach einem Gewicht durch Anwendung der

Kettenregel das in (7) angegebene δ ergibt.

Probleme und Lösungsansätze

Die Verwendung des Gradientenabstiegs bringt nicht nur Vorteile mit sich. Man muss sich

auch Gedanken um etwaige Probleme machen, die bei der Wahl dieser Methode auftreten


© modlab® 2007

können. Ohne Zweifel ist es mit Hilfe des Gradientenabstiegs möglich, das Netz immer besser

werden zu lassen. Es stellt sich nur die Frage, ob man nach Beendigung des Trainings

tatsächlich die optimale Gewichtskonstellation und dadurch auch die Funktion gefunden hat,

die die Datenklassen am besten voneinander trennen kann. Hierbei trifft man oft auf ein

allgemeines Problem von Optimierungsverfahren: das Auftreten lokaler Minima, die nicht

wieder verlassen werden.

Die Anzahl potentieller lokaler Minima erhöht sich mit der Dimension des Raumes. Es macht

also Sinn, das Netz nicht unbedingt mit der größtmöglichen Anzahl Parameter, die einem zur

Verfügung stehen, zu trainieren, sondern im Vorfeld wichtige Parameter zu selektieren.

Im Gegensatz dazu steht das Problem der flachen Plateaus. Diese bewirken, dass der Gradient

sehr niedrige Werte nahe 0 oder unter Umständen sogar 0 annimmt. Dies führt zu einer

Stagnation des Verfahrens. Das Risiko des Auftretens flacher Plateaus steigt mit abnehmender

Anzahl Trainingsparameter.

Im Wesentlichen sind zwei Methoden zu benennen, die sich mit den genannten Problemen

befassen, und gute Resultate erzielen:

1. Das Verändern der Lernrate zur Laufzeit, und

2. das Einführen eines Momentum-Terms.

Wie wir bereits gesehen haben, besitzt die Lernrate Einfluss auf die Veränderung der

Gewichte nach jedem Trainingszyklus. Am Anfang macht es durchaus Sinn, mit einer hohen

Lernrate zu starten, da der Algorithmus dadurch in der Lage ist, sich recht schnell auf dem

Suchraum zu bewegen. Im späteren Verlauf des Trainings ist eine große Lernrate allerdings

nicht von Vorteil. Diese stellt nämlich auch die Gefahr dar, gute Optima wieder zu verlassen

oder zwischen gleich guten Werten zu oszillieren. Eine Verringerung der Lernrate zur

Laufzeit des Algorithmus bietet den Vorteil, dass man zunächst recht schnell gute Positionen

des Raumes erkennt und diese durch die Verringerung der Lernrate nicht sofort wieder

verlässt. Weiterhin ist dadurch eine feinere Optimierung möglich, indem das Netz mit immer

kleiner werdenden Schritten weiter verbessert wird. Ein einfaches Beispiel für dieses

Vorgehen ist die Verwendung einer adaptiven Lernrate nach Darken und Moody. Diese ist

in (10) definiert.

r

nn

+

=

1

0ηη . (10)


© modlab® 2007

Man startet zunächst mit einer hohen Lernrate η0. Der Wert der Lernrate verringert sich im

Laufe des Trainings auf ungefähr η0 / (1 + n), wobei n die Anzahl bisher präsentierter

Trainingsbeispiele ist. Der Parameter r gibt hierbei die Stärke der Verringerung der Lernrate

in jedem Schritt an. Nach r Schritten hat sich die Lernrate halbiert. In der Praxis liefert bereits

die Verwendung dieses einfachen Verfahrens gute Ergebnisse.

Genauere Informationen zu dem Vorgehen nach Darken und Moody sowie eine Auswahl

weiterer möglicher Verfahren zur Veränderung der Lernrate während des Trainings sind in [4]

zu finden.

Die zweite erwähnte Möglichkeit, das Risiko lokaler Minima und flacher Plateaus zu

minimieren, ist die Einführung des sog. Momentum-Terms, der sich am ehesten durch

Impuls oder Trägheitsmoment übersetzen lässt. Dieses Konzept stammt aus der Physik. Im

Prinzip kann man sich den Suchraum (Fehlerfunktion) als zerklüftetes Gebirge vorstellen, mit

vielen Gipfeln und Tälern und möglicherweise auch Plateaus. Durch die Veränderung der

Gewichte bewegt man sich in diesem Suchraum von Gipfel zu Tal; der Fehler wird minimiert.

Durch die Stärke der Gewichtsänderung ist es möglich zu erkennen, wie stark das Gefälle ist.

Man stelle sich eine Kugel vor, die in dem vorhandenen Gebirge einen Hang hinunterrollt. Je

stärker das Gefälle ist, desto höher ist die Beschleunigung, die die Kugel erreicht. Erreicht die

Kugel das Tal, stoppt sie nicht sofort. Der Schwung bewirkt, dass die Bewegung erst

allmählich abklingt. Aufgrund dessen ist es mit Hilfe des Momentum-Terms möglich, flache

Plateaus zu überqueren oder zerklüftete Flächen zu passieren. Hierbei werden bei der

Berechnung der neuen Gewichtsveränderungen, die Änderungen des letzten Schrittes

berücksichtigt. Diese Änderungen wirken sich, mit einer gegebenen Stärke des Momentums,

auf die neuen Gewichtswerte aus. Mathematisch ist dies durch (11) definiert.

)1()( ,, −∆⋅+=∆ twtw jiijji αοδη , wobei (11)

)(, tw ji∆ die Gewichtsveränderung im aktuellen Schritt, )1(, −∆ tw ji die des vorherigen

Schrittes, jδ den Wert des unterhalb des Gewichtes liegenden Neurons j, iο die Ausgabe des

darüber liegenden Neurons i und α den Momentum-Parameter darstellen.

Wie man anhand von (11) erkennen kann, ist die Stärke des Momentums durch den Parameter

α modifizierbar. (11) stellt also eine um das Momentum erweiterte Form von (6) dar.

Trotz der offensichtlichen Vorteile, die der Momentum-Term mit sich bringt, sollte dennoch

darauf hingewiesen werden, dass dieser auch die Gefahr birgt, durch einen hohen Schwung


© modlab® 2007

aus einem besseren lokalen Optimum in ein schlechteres zu gelangen. Leider gibt es keinen

allgemeinen, problemunspezifischen Wert für Parameter α. Dessen Wert muss deswegen

durch trial-and-error Verfahren bestimmt werden.

Generell ist zu sagen, dass man beim Training neuronaler Netze eine Multistartstrategie

wählen und mehrere Netze mit verschiedenen Parametereinstellungen parallel trainieren

sollte. Die Selektion des besten Netzes ist dabei durch den RMSE (Abkürzung für root mean

square error) möglich. Dessen Berechnung ist in (12) gezeigt.

∑=

−=n

iii ot

nRMSE

1

2)(1

, wobei (12)

n die Anzahl der Trainingsbeispiele, ti der Zielwert und oi der berechnete Output-Wert des

Netzes für Trainingsbeispiel i ist.

Literatur

[1] W. Oberhofer, T. Zimerer 1996: „Wie Künstliche Neuronale Netze Lernen: Ein Blick in die Black Box der Backpropagation Netzwerke“ http://www.fh-ansbach.de/~thomas.zimmerer/pdfs/Literatur/DP_287.pdf [2] P. Baldi, S. Brunak 2001: „Bioinformatics - The Machine Learning Approach“, The MIT Press, Cambridge. [3] D. Kriesel 2005: „Ein kleiner Überblick über Neuronale Netze” http://www.dkriesel.com/content/studium/neuronetzefiles/neuronalenetze.pdf [4] W. Schiffmann, M.Joost, R.Werner: „Optimization of the Backpropagation Algorithm for Training Multilayer Perceptrons” http://www.cs.bham.ac.uk/~pxt/NC/schiffmann.bp.pdf

Partikelschwarm-Optimierung 63

© modlab® 2007

Partikelschwarm-Optimierung (PSO)

Die Idee für die PSO hat ihren Ursprung in der Simulation der Dynamik von biologischen

Systemen, wie Vogel- oder Fischschwärmen. Relativ früh schon wurde jedoch die Eignung

des Algorithmus als heuristisches Optimierungsverfahren erkannt und daraufhin

weiterentwickelt. Zum ersten Mal wurden Partikelschwärme von Eberhart und Kennedy im

Jahr 1995 als Optimierungsparadigma vorgeschlagen [1, 2].

Künstliche Schwärme bestehen aus Individuen, den so genannten Agenten oder Partikeln.

Indem den Partikeln einfache Regeln vorgegeben werden, nach denen sie sich zu verhalten

haben, können verschiedene Verhaltensmuster des Schwarmes entstehen. Trotz der

Einfachheit des Regelwerkes und somit der Einfachheit des Verhaltens der einzelnen

Individuen kann das Verhalten eines Schwarmes sehr komplex und vielschichtig sein. Solche

Komplexität entwickelt sich aus der Fähigkeit der Partikel, miteinander zu interagieren.

Indem sie untereinander Informationen über ihre Umwelt und über sich selbst austauschen,

beeinflussen sie sich auch gegenseitig in ihrem Verhalten.

Durch diesen Informationsaustausch können die Partikel gegenseitig ihre „Flugrichtung“ und

„Fluggeschwindigkeit“ im durch das Problem aufgespannten Suchraum (oder

Fitnesslandschaft) beeinflussen. Jeder Partikel hat das Bestreben, sich auf optimale, im Raum

gefundene Punkte hinzuzubewegen, so dass sich letztlich der Schwarm gerichtet durch den

Raum bewegt und nach einiger Zeit in einer Lösung konvergiert. Eine künstliche

Optimierungslandschaft ist in Abb. 1 gezeigt.


© modlab® 2007

Abbildung 1: Beispiel einer Fitnesslandschaft mit vielen lokalen Maxima und Minima (multimodale Fitnesslandschaft)

Der Algorithmus

Zu Anfang der PSO wird in dem von einem Problem aufgespannten Lösungsraum eine

bestimmte Anzahl von Partikeln generiert. Die Positionen dieser Partikel sind zunächst

zufällig gewählt, ebenso ihre Bewegungsrichtung und Geschwindigkeit. Jeder Partikel des

Schwarmes kann sich frei im Suchraum bewegen (Abb. 2).

Abbildung 2: Momentaufnahme: Partikel in einer Fitnesslandschaft, das globale Optimum ist durch den Schnittpunkt der schwarzen Geraden mit der Lösungsoberfläche visualisiert.

Bei einem D-dimensionalen Problem wird die Position jedes Partikels – hier des i-ten

Partikels – durch einen Ortsvektor Xi = (xi1, xi2, xi3,…xiD) beschrieben. Die Geschwindigkeit

des i-ten Partikels wird analog mit einem Geschwindigkeitsvektor Vi = (vi1, vi2, vi3,…viD)

beschrieben. Da jedem Partikel durch den Ortsvektor eine Position im Lösungsraum

Dimension 1

Funktio

nsw

ert

Dim

ension 2


© modlab® 2007

zugewiesen ist, ist es auch möglich, die jeweilige Fitness der Partikel in Bezug auf das

Problem zu bestimmen. Des Weiteren existieren zwei Arten von „Gedächtnis“ in der PSO:

• Zum einen haben Partikel ein so genanntes „individuelles“ Gedächtnis, in dem sie

die beste Position pi abspeichern können, an der sie sich jeweils selbst im

Optimierungsprozess bis zum aktuellen Zeitpunkt befunden hat.

• Zum anderen besitzen die Partikel ein „soziales“ Gedächtnis, das man sich auch als

„Schwarmgedächtnis“ vorstellen kann. Dieses „soziale“ Gedächtnis speichert die

Position pb ab, die bis zum aktuellen Zeitpunkt die insgesamt beste, von allen Partikeln

zusammen gefundene Position ist.

Ausgehend von diesen beiden Informationen – also den Positionen aus den individuellen

Gedächtnissen und dem sozialen Gedächtnis – sowie dem aktuellen Geschwindigkeitsvektor

werden nun für jeden Partikel neue Geschwindigkeitsvektoren berechnet.

In dem 1998 von Shi und Eberhart vorgeschlagenen PSO-Algorithmus [3] erfolgt die

Berechnung nach Gleichung 1:

))(())(()(v)1(v 2211ii txprntxprntwt ibii −⋅⋅+−⋅⋅+⋅=+ , (1)

wobei t für die Epochenzahl steht, n1 und n2 Konstanten zur Gewichtung des Einflusses der

beiden Gedächtnisse auf das Schwarmverhalten sind, und r1 und r2 Zufallszahlen zwischen 0

und 1 sind. Das so genannte „Massenträgheitsgewicht“ w gewichtet den Einfluss des alten

Geschwindigkeitsvektors vi eines Partikels und dient maßgeblich der Kontrolle der

Gesamtgeschwindigkeit, mit der sich Partikel durch den Suchraum bewegen.

„Vernünftige“ Werte für die Konstanten, die jedoch von Anwendung zu Anwendung

unterschiedlich passend sein können, liegen für n1 und n2 zwischen 1 und 2, für w zwischen

0,5 und 1. Diese Angaben sind jedoch nur als Anhaltspunkte zu verstehen.

Um einer unkontrollierten, explosionsartigen Bewegung des Schwarmes vorzubeugen, wird

oftmals der maximale Wert, den der Geschwindigkeitsvektor annehmen kann, eingeschränkt.

Wird dieser maximale Wert überschritten, so wird vi auf den Grenzwert zurückgesetzt. Eine

weitere Möglichkeit, die Geschwindigkeit des Schwarmes zu kontrollieren, ist das

Massenträgheitsgewicht mit zunehmender Epochenzahl linear kleiner werden zu lassen, was

in der Regel zu einer Konvergenz des Schwarmes über die Zeit führt [4].


© modlab® 2007

Die Positionen der einzelnen Partikel werden durch Addition des in jeder Epoche neu

errechneten Geschwindigkeitsvektors zu dem Ortsvektor neu berechnet (Gleichung 2):

)1()()1( ++=+ tvtxtx iii . (2)

Nach dem Update der Positionen kann ein neuer Evaluationszyklus beginnen, bis ein

Stoppkriterium (wie z.B. eine maximale Anzahl an Epochen) erreicht wird.

In Pseudocode lässt sich die PSO folgendermaßen beschreiben:

begin

Positionen und Geschwindigkeiten der Partikel initialisieren;

Fitness evaluieren;

Gedächtnisse initialisieren;

while (Abbruchbedingung ≠ true)

neue Geschwindigkeiten berechnen;

Positionen neu berechnen;

Fitness evaluieren;

Gedächtnisse auffrischen;

end

end

Zusammenfassung

Die Partikelschwarm-Optimierung ist ein biologisch motiviertes, leicht zu implementierendes

und effizientes heuristisches Optimierungsverfahren. Im Vergleich zu anderen biologisch

motivierten Optimierungsalgorithmen, wie beispielsweise genetischen Algorithmen [5], hat es

nur wenige freie Parameter, die jedoch entscheidenden Einfluss auf das Verhalten der

Schwärme haben.

PSO wurde in der Vergangenheit bereits erfolgreich für unterschiedlichste Anwendungen

eingesetzt, vom Training künstlicher neuronaler Netze bis hin zum Lösen technischer

Probleme im Ingenieurswesen.


© modlab® 2007

Literatur

1. Kennedy J, Eberhart RC: Particle swarm optimization. In Proceedings of IEEE International Conference on Neural Networks; Piscataway, NJ. 1995: 1942-1948.

2. Eberhart RC, Kennedy J: A new optimizer using particle swarm theory. In

Proceedings of the Sixth International Symposium on Micromachine and Human Science; Nagoya, Japan. 1995: 39-43.

3. Shi Y, Eberhart RC: Parameter selection in particle swarm optimization. In

Evolutionary Programming VII: Proceedings of the Seventh Annual Conference on Evolutionary Programming; New York, USA. 1998: 591-600.

4. Kennedy J, Eberhart RC: Swarm Intelligence. San Diego: Academic Press; 2001. 5. Holland J.H.: Adaptation in natural and artificial system, Ann Arbor, The University

of Michigan Press; 1975.

Ameisenalgorithmen 69

© modlab® 2007

Ameisenalgorithmen

Ameisenalgorithmen stammen aus dem Gebiet der naturanalogen bzw. biologisch motivierten

Algorithmen. Das Ant Colony Optimization (ACO) Konzept wurde von Marco Dorigo

(Dorgio et al., 1996) als Metaheuristik eingeführt (http://www.aco-metaheuristic.org/). Die

Idee des Ameisenalgorithmus wurde dem Foraging (Futtersuche) bestimmter

Ameisenunterfamilien nachempfunden. Diese sind in der Lage, bestehende Hindernisse

zwischen Nest und Futterquelle bei verschiedenen möglichen Wegen kurzmöglichst zu

verbinden (Deneubourg et al., 1990; Abb.1).

Abbildung 1. Zeitlicher Ablauf der Wegfindung bei Ameisen. Die Graustufe der Pfeile spiegelt die Pheromonkonzentration wieder.

Die Ameisen verwenden zuerst alle zur Verfügung stehenden Verbindungen zwischen

Futterquelle und Nest (Abb. 1a). Über die Zeit (Abb. 1b) ist ein autokatalytischer Prozess zu

beobachten, der immer mehr Arbeiterinnen auf den kürzeren Weg zieht, bis schlussendlich

(Abb. 1c) alle Wege außer dem kürzesten ignoriert werden. Diese Lösung eines natürlichen

Optimierungsproblems ist unter den folgenden Gesichtspunkten für die Informatik interessant:


© modlab® 2007

• Keine übergeordnete Intelligenz (Königin) steuert die Ameisen.

• Einzelne Ameisen haben kein Problembewusstsein.

• Die optimale Lösung wird gefunden.

• Die Ameisen kooperieren, um eine Lösung zu finden.

Der Aspekt der Kooperation ist von entscheidender Bedeutung. Den Ameisen ist es möglich,

über abgegebene Duftstoffe (Pheromone) über ihre Umwelt als Medium zu kommunizieren.

Diese Kommunikation über eine Veränderung der Umwelt, um einem gemeinsamen Ziel zu

folgen, bezeichnet man als Stigmergie, eingeführt von Grassé 1959 für Termiten (Grassé, P.

P. & Noirot, C., 1959; Theraulaz, G. & Bonabeau, E., 1999). Die Ameisen lösen das in Abb.1

gezeigte Weg-Problem mit Hilfe von Pheromonen (Deneubourg et al., 1990). Auf dem

kürzeren Weg wandern pro Zeit mehr Ameisen als auf dem längeren Weg, wodurch die

Pheromonkonzentration auf diesem schneller anwächst (bzw. weniger abnimmt). Die aus dem

Nest nachfolgenden Ameisen werden durch die Pheromonkonzentration hinsichtlich ihrer

Entscheidung für einen Weg zur Futterquelle beeinflusst. Die kürzere Route wird durch die

höhere Pheromonkonzentration mit einer gewissen Wahrscheinlichkeit bevorzugt.

Einzelne gleichartige Individuen werden also über eine Kommunikationsschnittstelle

verbunden und dadurch in die Lage versetzt, Probleme zu lösen. Es entsteht eine Lösung aus

dem System, die über die Fähigkeiten der einzelnen Ameisen hinausgeht (Emergenz). Dies

wird durch die Ameisenalgorithmen für die Informatik im Bereich von

Optimierungsproblemen adaptiert. Ameisen werden formal zu „Agenten“. Der Aspekt der

Pheromone wird 1:1 übernommen.

Ameisenalgorithmen in der Anwendung

ACS (Ant Colony System) für das Travelling Salesman Problem nach Dorigo &

Gambardella (1997).

Im Problem des Handlungsreisenden (TSP) geht es darum, eine möglichst kurze Rundreise

zwischen verschiedenen Städte zu organisieren, wobei jede Stadt nur einmal besucht werden

darf. Ein Beispiel ist in Abb.2 zu sehen.


© modlab® 2007

Abbildung 2. Vollständiger verbundener Graph für das TSP Problem. Alle Kanten haben unterschiedliche Gewichte (zufällig) alle Knoten eindeutige Bezeichnungen.

In Abb.2 stellen Knoten, Städte und Kanten, Straßen mit unterschiedlichem Gewicht (Länge)

entsprechen. Ein Agent (k) startet nun in einer zufälligen Stadt und wählt von dieser Stadt (r)

die nächste noch nicht besuchte Stadt (s) entsprechend Gleichung (1).

( )[ ] ( )[ ]{ }

≤⋅= ∉

sonstS

qqwennururs kMu

0,,maxarg βητ

, wobei (1)

( )ur,τ der Pheromonkonzentration auf der Kante zwischen Stadt r und u entspricht. ( )ur,η

ist eine heuristische Funktion, invers zu der Distanz zwischen den Städten r und u. β

gewichtet die relative Bedeutung der Pheromonspur sowie die Entfernung (Kantengewicht)

zwischen den Städten. q wird mit gleichverteilter Wahrscheinlichkeit aus [0,1] gezogen, 0q

ist ein Parameter, S ist eine zufällig nach Gleichung 2 belegte Variable.

( )

( )[ ] ( )[ ]( )[ ] ( )[ ]

∉

⋅

⋅

= ∑∉

sonst

Mswennurur

srsr

srpk

Muk

k

0

,,

,,

,β

β

ητ

ητ

. (2)

Die Formel bevorzugt Kanten mit einer stärkeren Pheromonspur und geringerem Gewicht

(Länge). Die Idee ist, dass der Agent mit einer Wahrscheinlichkeit von 0q die Möglichkeit

hat, auf das kollektive Gedächtnis zuzugreifen, oder mit einer Wahrscheinlichkeit von

01 q− eine „eigene“ Stadt nach den Vorgaben von Gleichung 2 zu wählen. Für eine

weiterführende Beschreibung des Pheromonupdates (globales, lokales) sei auf Dorigo &

Gambardella (1997) verwiesen.


© modlab® 2007

Ein Applet für TSP nach dem ACO Konzept ist hier zu finden:

http://www.geocities.com/andres_perezuribe/

Es ist zu beachten, dass das hier verwendete Pheromonupdate explizit auf die Problemstellung

zugeschnitten ist. Das ACO Konzept lässt sich jedoch auf eine weites Spektrum von

Optimierungsproblemen anwenden. Für eine Auflistung möglicher Anwendungen sei auf

Dorigo & Stützle, 2004 verwiesen. Die Herausforderung bei der Anwendung von ACO liegt

darin,

a) eine geeignete Methode des Pheromonupdates zu finden, und

b) das „Pfad“-Konzept an die Problemstellung anzupassen.

Naheliegende Anwendungen sind Verkehrssimulationen oder Routing Probleme. Auch die

Möglichkeit der Modellierung biologischer Systeme (echter Ameisen) stellt eine mögliche

Anwendung da.

Literatur

Deneubourg, J-L., Aron, S., Goss, S. & Pasteels, J.M. (1990). The self-organizing exploratory pattern of the Argentine ant. Journal of Insect Behavior 3: 159-168. Dorigo, M., Maniezzo, V. & Colorni, A. (1996). Ant system: optimization by a colony of cooperating agents. IEEE Transactions on Systems, Man, und Cybernetics-Part B. 26: 29-41. Dorigo M. & Gambardella, L.M. (1997). Ant Colonies for the Traveling Salesman Problem. BioSystems 43:73-81. Also Tecnical Report TR/IRIDIA/1996-3, IRIDIA, Université Libre de Bruxelles. Dorigo, M. & Stützle, T. (2004). Ant Colony Optimization. The MIT Press, Cambridge, Massachusetts. Grassé, P. P. & Noirot, C. (1959). The development of symbiosis in Isoptera. Experientia 15: 365-72. Theraulaz, G. & Bonabeau, E. (1999). A brief history of stigmergy. Artificial Life 5: 97-116.

lab Grundlagen der Bioinformatik mod - Benutzer-Homepageskrausk/Bioinformatik/Bioinfo... · BLAST...

Documents

Transcript of lab Grundlagen der Bioinformatik mod - Benutzer-Homepageskrausk/Bioinformatik/Bioinfo... · BLAST...