Phylogenetische Netzwerke

 

Dissertation

zur Erlangung des Doktorgrades

des Fachbereichs Mathematik

der Universität Hamburg

 

 

 

Vorgelegt von

Arne Röhl

aus Pinneberg

 

 

Hamburg

1999

 

 

 

 

 

Als Dissertation angenommen vom Fachbereich

Mathematik der Universität Hamburg

  

auf Grund der Gutachten von Prof. Dr. H. J. Bandelt

und Prof. Dr. B. Brinkmann

  

Hamburg, den 03.06.1999 

  Prof. Dr. H. Daduna

Dekan des Fachbereichs Mathematik

 

 

 

 

Inhaltsverzeichnis

 

Einleitung *

Phylogenetische Netzwerke *

Sparsamste (MP) Bäume / Steinerbäume *

Medianvektoren und Quasimediane *

Quasimediane Erzeugung *

Steinerbäume/ sparsamste (MP) Bäume in quasimedianen Netzwerken *

Netzwerk- Algorithmen *

Median Joining Algorithmus (MJ) *

Reduced Median Algorithmus (RM) *

Vergleich von MJ mit RM *

Das Verfahren von Foulds et al. (1979) / greedy FHP *

Vergleich von MJ mit dem greedy FHP-Verfahren *

Vergleich von RM mit dem greedy FHP-Verfahren *

Wahl eines Verfahrens *

Star Contraction (SC) Algorithmus *

Algorithmus *

Vergleich mit Wills (1995) *

Beispiel: Chen et al. (1995) *

Effektivitäts- und Laufzeitanalyse *

Hochdimensionale Netzwerke / Cycle Pruning (CP) Algorithmus *

CP- Algorithmus *

Verfahren von Excoffier et al. (1994) *

Vergleich des CP-Algorithmus mit dem Verfahren von Excoffier et al. (1994) *

Beispiel: Excoffier et al. (1994) *

Simulationsstudie *

Vorgehensweise *

Datensatz 1: *

Datensatz 2: *

Anhang A: NETWORK 2.0 *

Datenformate *

MJ- Algorithmus (Bandelt et al. (1999)) *

RM-Algorithmus (Bandelt et al. (1995)) *

Graphikprogramm *

Star Contraction (SC) Algorithmus *

Cycle Pruning (CP) Algorithmus *

Literatur *

Anhang B: Rolf B, Röhl A, Forster P, Brinkmann B. On the Genetic origins of the Turks *

Zusammenfassung *

Lebenslauf *

 

Einleitung

 

Der Ursprung und die Entwicklung der Menschheit ist seit Jahren ein vieldiskutiertes Thema. Die klassische Anthropologie hat sich zur Erforschung dieses Gebietes auf äußerliche Merkmale wie Hautfarbe, Körpergröße und ähnliches beschränkt und mit der geographischen Nähe der Verbreitungsgebiete kombiniert. Hierbei treten allerdings zwei Probleme auf:

 Ausweg aus der Problematik, ob äußerliche Unterschiede aufgrund der Vererbung oder aber der Umwelteinflüsse bestehen, bietet ein neuerer Ansatz, der zusätzlich genetisches Material zur Postulierung von Abstammungsverhältnissen benutzt: Die Erbsubstanz (DNS), bestehend aus den vier Nukleinbasen Adenin (A), Guanin (G), Thymin (T) und Cythosin (C), angeordnet auf einer Doppelhelix, die die Merkmale und Eigenschaften eines jeden Lebewesens bestimmt. Äußerliche Anpassungen der Menschen an Umwelteinflüsse haben (soweit bekannt) nur Einfluß auf einen geringen Teil der DNS, so daß aus Ähnlichkeiten der DNS auch auf Ähnlichkeiten der Populationen geschlossen werden kann.

Der weitaus größte Teil der Erbsubstanz wird dabei in jeder Generation aus väterlicher und mütterlicher DNS rekombiniert und eignet sich deshalb insbesondere dazu, Verwandtschaftsverhältnisse innerhalb von Familien festzustellen und in der Forensik anhand des ‘genetischen Fingerabdruckes’ Verbrecher zu überführen (Jeffreys et al. (1985)). Die Rekonstruktion einer Abstammungsfolge und damit einer Zeitskala ist aber mit dieser Art der DNS nur für kleine Genabschnitte möglich, die normalerweise zu kurz sind, um detaillierte Aussagen bzgl. der Ausbreitung des modernen Menschen machen zu können. Bestenfalls können Wahrscheinlichkeiten für den geographischen Ursprung von Individuen berechnet werden (Meyer et al. (1995)).

Es gibt allerdings auch zwei nicht rekombinierende Bereiche der DNS, die aufgrund dieser Eigenschaft bzgl. der phylogenetischen Analyse von Populationen geeignet sind: Zum einen die DNS der Mitochondrien (mtDNS), ein kleines, im Menschen 16.569 Basenpaare (bp) langes Ringmolekül, das für einige Stoffwechselvorgänge zuständig ist, zum anderen das Y-Chromosom. Diese DNS wird ausschließlich von einer elterlichen Seite vererbt (die mtDNS mütterlicherseits und das Y-Chromosom väterlicherseits), so daß eine Neukombinierung der DNS verschiedener Individuen nicht möglich ist.

Allein diese Tatsache wäre aber noch nicht ausschlaggebend für die Eignung dieser speziellen DNS als genetisches Material zur Bestimmung von Abstammungs-verhältnissen, da die Menschheit aus einem Genpool hervorgegangen ist, und insofern keinerlei Unterschied zwischen verschiedenen Individuen festzustellen wäre, gäbe es nicht die Eigenschaft der DNS-Polymerase, bei der Kopie der DNS kleine Fehler (Mutationen) einzubauen. Zu unterscheiden sind dabei Transversionen (A oder G mutiert nach C oder T bzw. umgekehrt) und Transitionen (C mutiert nach T oder A nach G bzw. umgekehrt), da Transversionen in der mtDNS seltener auftreten als Transitionen. Einmal aufgetretene Mutationen werden dann von Generation zu Generation weitervererbt, so daß es theoretisch möglich ist, den exakten Stammbaum der aufgetretenen Mutationen zu rekonstruieren. Da außerdem die Mutationen einigermaßen regelmäßig auftreten, ist auch der Aufbau einer Zeitskala realisierbar, anhand derer Aussagen über Zeitpunkte möglich werden. Im Bereich der Nukleotidpositionen (np) 16090- 16369 des hypervariablen Segments I der 1127 bp langen Kontrollregion der mtDNS, die verglichen mit der kompletten mtDNS hochmutabel ist (Horai et al. 1995), rechnet man ungefähr mit einer Transition alle 20.180 Jahre (Forster et al. (1996)).

Im Y-Chromosom muß zusätzlich zwischen Punktmutationen und Short Tandem Repeats (STRs) unterschieden werden. Bei letzteren Merkmalen ändert sich nämlich nicht nur wie bei Punktmutationen ein Nukleotid, sondern eine kurze spezifische Sequenz von Nukleotiden tritt in unterschiedlicher Anzahl hintereinander auf. Verglichen mit Punktmutationen ist die Mutationsrate bzgl. der Anzahl der Repeats wahrscheinlich recht hoch (Weber et al. (1993)), so daß diese Art von Merkmalen sich insbesondere zur Untersuchung nah miteinander verwandter Populationen eignet (Rolf et al., vgl. Anhang B).

Den Vorgang einer Mutation kann man sich so vorstellen wie ein Buch, welches immer wieder abgeschrieben wird. Ein einmal vorhandener Fehler wird weiterkopiert und evtl. werden Fehler hinzugefügt, so daß letztlich ein klarer Baum der Abstammungshierarchie aufgestellt werden kann. Würde man außerdem noch das Originalwerk kennen (die Wurzel), so könnte man einen Stammbaum aufstellen. Ein Problem gibt es jedoch: Sollte bei einer Kopie ein Fehler rückgängig gemacht werden bzw. unabhängig voneinander der gleiche Fehler bei verschiedenen Kopien auftreten, weisen Kopien Ähnlichkeiten zueinander auf, die aufgrund der wahren Abstammungsfolge eigentlich gar nicht vorhanden sind (Homoplasien). Eben diese Problematik tritt auch in der mtDNS und bei Y-STRs auf. Rückmutationen können auftreten und aufgrund der relativ kurzen Sequenz der mtDNS werden insbesondere in den hypervariablen Segmenten I und II der Kontrollregion (HVS I, HVS II) auch unabhängig voneinander Mutationen an der gleichen Stelle auftreten. Bei den Y-STRs kommt die Problematik hinzu, daß bisher ungeklärt ist, ob die Anzahl der Repeats an einem Locus nach oben beschränkt ist, so daß bei Erreichen der Grenze nur noch Rückmutationen auftreten können. Außerdem kann sich die Länge eines STRs bei einer Mutation gleich um mehrere Repeats verlängern oder verkürzen und es ist ungewiß, wie oft dies auftritt.

 

Phylogenetische Netzwerke

 

Um Homoplasien gerecht zu werden, gibt es prinzipiell zwei Möglichkeiten: zum einen können bei unklaren Abstammungsverhältnissen mehrere plausibel erscheinende Bäume als Lösungen in Betracht gezogen werden, zum anderen besteht die Möglichkeit, sämtliche Lösungen (plausible Bäume) in einem Netzwerk, d. h. einem Lösungsraum, zu vereinigen. Liegen beispielsweise nach wiederholter Kopie eines Textes und der Kopien des Textes die drei fehlerhaften Worte ‘Leerer’, ‘Lehre’ und ‘Leere’ des ursprünglichen Wortes ‘Lehrer’ vor und legt man zugrunde, daß nicht gleichzeitig zwei Fehler aufgetreten sind (‘one-step mutation model‘), so können die Worte auf mehrere Arten aus dem ursprünglichen Wort entstanden sein:

 

 

Noch schwieriger wird das Problem, wenn nur die vier Worte vorliegen, aber nicht bekannt ist, welches das ursprüngliche (die Wurzel) ist. Insofern ist zu unterscheiden zwischen gerichteten und ungerichteten Bäumen. Bei letzteren ist nicht bekannt, welchen Ursprung der Baum hat, so daß keine Mutationsrichtungen angegeben werden können.

Da die genaue Abstammungsfolge der Worte voneinander unbekannt ist (bekannt ist in diesem Beispiel nur, daß ‚Lehrer‘ das Originalwort ist), ist es zweckmäßiger und übersichtlicher, alle vier bzgl. der Annahme des ‚one-step mutation model‘ gleich guten Lösungen in einem Netzwerk zu vereinigen.

Nach Festlegung des Wurzelknotens kann jetzt durch Entfernung einer Kante jeder der zuvor gezeigten Bäume generiert werden.

Aber auch diese Art der Darstellung kann problematisch sein, falls in einem Datensatz sehr viele dieser Inkompatibilitäten von Merkmalen (in obigem Beispiel die Buchstaben 3 und 6) vorliegen, insbesondere wenn mehrere Merkmale paarweise inkompatibel sind. Beispielsweise kann bei drei paarweise inkompatiblen Merkmalen ein dreidimensionaler und bei vier derartigen Merkmalen ein vierdimensionaler Würfel (vgl. Abb. 1) entstehen. Das resultierende Netzwerk kann schlechtestenfalls zu einem hoch-dimensionalem Gebilde ‘entarten’, das kaum zeichen- und deutbar ist. Praktische Erfahrungen haben aber gezeigt, daß dies bei Datensätzen der mtDNS selten auftritt (Bandelt et al. (1995), Bandelt et al. (1999), Watson et al. (1997), Forster et al. (1996)), und sich diese Ereignisse in der Regel durch fehlerhafte Daten und hochmutable Positionen erklären lassen, so daß normalerweise die Netzwerkpräsentation bei intraspezifischen Daten der Baumdarstellung vorzuziehen ist, zumal ein Baum jederzeit aus dem Netzwerk durch Zusatzinformationen bzw. geeignete Verfahren herausgefiltert werden kann, indem aus jedem im Netzwerk enthaltenen Kreis genau eine Kante entfernt wird.

 

 

Fig. 1: Vierdimensionaler Hyperwürfel, aufgespannt durch die fünf schwarz gekennzeichneten Eckpunkte.

 

Zu bedenken ist außerdem, daß bei Anwendung eines Baumverfahrens wie Neighbour Joining (NJ) oder Maximum Likelihood (ML) auf einen Datensatz, der Strukturen unterstützt, wie die in Fig. 1 gezeigte, nur ein Bruchteil der enorm vielen optimalen Lösungen als Ergebnis ausgegeben wird.

 

Sparsamste (MP) Bäume / Steinerbäume

 

Zusätzliche Probleme bei der Rekonstruktion können auftreten, falls, um auf das Beispiel des Buches zurückzukommen, einige Kopien des Buches fehlen, entweder weil diese Kopien nicht mehr existieren, oder weil nicht alle Kopien in der Sammlung enthalten sind. Wurden dann von einem dieser fehlenden Bücher mehrere fehlerhafte Kopien erstellt, so ist in dem entsprechenden Baum der zugehörige Verzweigungspunkt nicht enthalten.

Unter der Annahme, daß bei der Kopie der gleichen Vorlage (bzw. der DNS) selten der gleiche Fehler unabhängig voneinander gemacht wird, sollten solche wahrscheinlich fehlenden Daten (sei es ein Buch oder eine DNS-Sequenz) ergänzt werden. Diese Ergänzung fehlender Verzweigungspunkte des Sequenzraumes ist in der Biologie bekannt als die Suche nach kürzesten Bäumen (MP-Bäumen) und mit dem ‘Steiner-Problem’ der Mathematik in einem metrischen Raum gleichzusetzen. Im wesentlichen stellt sich die Problematik so dar: Gegeben ist eine Anzahl Knoten (Bücher, DNS-Sequenzen,...) und Distanzen der Knoten zueinander (z. B. die Hamming-Distanz, die mit der Anzahl der Unterschiede gleichzusetzen ist), so daß ein metrischer Raum vorliegt. Durch Hinzunahme geeigneter Zwischenknoten soll ein Baum minimaler Länge (Steinerbaum) bzgl. der gegebenen Distanzen generiert werden.

 

Medianvektoren und Quasimediane

 

Für drei gegebene Sequenzen ist es möglich, spezielle Lösungen dieses Problems durch die quasimediane Operation (vgl. Bandelt et al. (1994)) zu finden. Ein Quasimedian ist eine Menge von Knoten (Medianvektoren) {u=(xyz), v=(yzx), w=(zxy)}, gebildet aus drei Ausgangssequenzen x, y, z Î X, mit folgenden Eigenschaften:

 

  1. u, v und w liegen auf kürzesten Wegen zwischen x, y und z und
  2. d(u, v) = d(u, w) = d(v, w) = k, wobei der Wert k so klein wie möglich sein soll.

 

d(u,v) bezeichne dabei die Distanz zweier Knoten u und v, d. h. die Länge eines kürzesten Weges zwischen u und v.

Besteht ein Quasimedian nur aus einem Medianvektor, entsprechend k=0, so heißt dieser Median des Tripels. Ein Medianvektor u=(xyz) nimmt dabei, bezogen auf Sequenzen der mtDNS, für jedes Merkmal den Majoritätszustand an. Ist dieser nicht eindeutig, d. h. alle drei Sequenzen weisen unterschiedliche Ausprägungen des Merkmals auf, so wird der Merkmalszustand der ersten Sequenz übernommen, in obigem Beispiel also der Merkmalszustand des Merkmales x.

Das allgemeine Steinerbaum-Problem gehört allerdings zur Klasse der NP-vollständigen Probleme (vgl. Karp (1972), Gusfield (1997)), d. h. es wurde bisher kein Algorithmus gefunden (und es spricht einiges dafür, daß ein solcher Algorithmus auch nicht existiert), der bei beliebiger Anzahl Ausgangssequenzen garantiert eine Lösung in polynomialer Zeit findet. Lawler konnte 1976 immerhin einen Algorithmus angeben, der in der Anzahl der Sequenzen polynomial anwächst und nur in der Anzahl der Steinerpunkte exponentiell. Daher müssen heuristische Verfahren benutzt werden, die zwar nicht mehr garantieren können, eine optimale Lösung zu finden, aber im allgemeinen dennoch dem Optimum recht nahe kommen.

 

Quasimediane Erzeugung

 

Gegeben sei ein endliches Alphabet L (im Falle von DNS Sequenzen gleichzusetzen mit L ={A, G, C, T}) und eine Menge von Sequenzen X Í L n, wobei n die Sequenzlänge bezeichne. Mit xi sei für alle i mit 1 £ i £ n die i-te Koordinate der Sequenz x Î X bezeichnet. In der gleichen Weise sei xJ die Beschränkung der Sequenz x auf alle Koordinaten der Menge Æ ¹ J Í {1,...,n}. Für eine Menge X setze man XJ:={xJ | x Î X}.

Da X eine endliche Menge ist, gibt es einen Index k ³ 0, für den mit der Abkürzung X[ 0 ]:=X und X[ i ]:={(xyz)| x, y, z Î X[ i - 1 ]} , i ³ 1, gilt: Q=Xk], wobei Q die von X erzeugte quasimediane Algebra ist.

Lemma 1: Sei Æ ¹ X Í L n, n ³ 2. Genau dann ist Q die von X erzeugte quasimediane Algebra, wenn QJ für alle zwei-elementigen Teilmengen J der Menge {1,...,n} die durch XJ erzeugte quasimediane Algebra in L 2 ist.

Beweis:

 

‘=>‘: Da die quasimediane Operation koordinatenweise vorgeht, muß jede Projektion einer Sequenz aus der durch X erzeugten quasimedianen Algebra Q auf die Koordinaten einer Teilmenge J der Koordinatenmenge {1,...,n} automatisch auch in der durch XJ erzeugten quasimedianen Algebra QJ liegen, folglich auch für jede zwei-elementige Teilmenge.

 

‘<=‘: Beweis durch Induktion nach der Anzahl der Elemente #J aus J:

QJ sei für alle zweielementigen Teilmengen J der Menge {1,...,n} die durch XJ erzeugte quasimediane Algebra in L 2. Nun gelte #J³3 und jede Beschränkung der Ausgangssequenzen aus X auf eine (#J – 1) - elementige Koordinatenteilmenge erzeuge die entsprechende Projektion der durch X erzeugten quasimedianen Algebra. Für eine Sequenz v Î Q seien i, j, k beliebige voneinander verschiedene Koordinaten aus J Í {1,...,n}. Dann gibt es Sequenzen w, x, y aus Q mit den folgenden Eigenschaften:

 

vJ - { i } = wJ - { i } ; vJ - { j } = xJ - { j } ; vJ - { k } = yJ - { k }

 

Insbesondere gilt dann:

 

vi = xi = yi ; vj = wj = yj ; vk = wk = xk und

vl = wl = xl =yl für alle lÎ J - {i, j, k}

 

Somit ist vJ=(wxy) J und wird durch X J erzeugt. Dies beendet den Induktionsschluß.

Q.e.d.

 

Um festzustellen, ob eine Sequenz in der durch vorgegebene Sequenzen X erzeugten quasimedianen Algebra Q liegt, reicht es also, dies für alle zwei-elementigen Teilmengen der Koordinaten zu überprüfen.

 

Lemma 2: Sei Æ ¹ X Í L 2. Eine Sequenz v Î L 2 gehört genau dann zu der durch X erzeugten quasimedianen Algebra Q, wenn eine der folgenden Bedingungen erfüllt ist mit {i , j } = {1, 2}:

 

(0) v Î X ;

(1) es gibt Sequenzen w, x, y in X mit vi = wi, vj = xj = yj und xi ¹ yi ;

(2) es gibt Sequenzen w, x, y, z in X mit vi = wi, vj = xj , yi = zi , xi ¹ yi und yj ¹ zj ;

(3) es gibt Sequenzen w, x, y, z in X mit vi = wi, vj = xj , wj = zj , xi ¹ yi und yj ¹ zj .

 

Die Sequenz v ist im Fall (k) in X[ k ] enthalten und es werden genau k quasimediane Operationen benötigt, um die Sequenz v zu generieren ( k Î {0,1, 2, 3} ).

 

Beweis: ‘=>‘: Sei v eine Sequenz der durch X erzeugten quasimedianen Algebra Q:

Angenommen, der Fall (0) tritt nicht ein, d. h. die Sequenz v ist nicht in X enthalten. Da das kartesische Produkt X1 x X2 der Projektionen von X auf die erste und zweite Koordinate eine quasimediane Algebra bildet, gilt Q Í X1 x X2. Somit existieren Sequenzen w und x in X mit v1 = w1 und v2 = x2 und notwendigerweise w1 ¹ x1 und w2 ¹ x2. Gibt es eine weitere Sequenz yÎX - {w, x}, die mit v in einer Koordinate übereinstimmt, z. B. v2 = y2 , dann gilt (1) mit i = 1 und j = 2. Anderenfalls (mit v1 = y1) gilt (1) mit i = 2 und j = 1 bei Vertauschen der Rollen von w und x.

Sei nun angenommen, w und x wären die einzigen Sequenzen in X, die in einer Koordinate mit v übereinstimmen. Hätten alle Sequenzen in X verschiedene Ausprägungen der ersten Koordinate sowie unterschiedliche Ausprägungen der zweiten Koordinate, so wären Q und X im Widerspruch zur Annahme identisch, daß v nicht in X enthalten ist. Demnach gibt es ein Paar von Sequenzen y¹ z in X – {w, x}, die an einer Koordinate gleiche Ausprägungen aufweisen, z. B. y1 = z1. Gilt außerdem x1 ¹ y1, so ergibt sich (2). Gilt dagegen y2 = z2 und w2 ¹ y2, dann gilt (2) nach Vertauschen der Rollen von w und x.

Im letzten Fall gibt es eine Sequenz z aus X – {w, x}, die in einer Koordinate mit w oder x übereinstimmt, beispielsweise z2 = w2. Bestünde X nur aus den Sequenzen w, x und einer Teilmenge der quasimedianen Subalgebra ({x1} x L ) È (L x {w2}), dann könnte v nicht in der durch X erzeugten quasimedianen Algebra enthalten sein, denn die Menge {{x1} x L ),(L x {w2}} ist bezüglich der quasimedianen Operation abgeschlossen. Somit gibt es eine Sequenz y aus X außerhalb der quasimedianen Subalgebra und (3) ist erfüllt.

 

‘<=‘: Um die Rückrichtung zu zeigen, braucht man nur nachzurechnen, daß gilt

v = (wxy) wenn (1) erfüllt ist,

v = (wx(xyz)) wenn (2) erfüllt ist und

v = (wx(xy(yzw))) wenn (3) gilt. Q.e.d.

 

Zur Überprüfung der Zugehörigkeit einer Sequenz zu einer durch eine Menge X Í L 2 generierten quasimedianen Algebra Q kann nun der folgende Algorithmus benutzt werden:

Man zähle für alle Ausprägungen der ersten und zweiten Koordinate die Anzahl der im Datensatz X vorkommenden Sequenzen, die die entsprechende Ausprägung an der jeweiligen Koordinate aufweisen.

Gibt es entweder keine Sequenz mit der Ausprägung v1 in der ersten Koordinate oder mit der Ausprägung v2 in der zweiten Koordinate, so folgt vÏ Q. Ist dies nicht der Fall, und die Summe der Sequenzen xÎ X mit Ausprägung x1 = v1 oder x2 = v2 ist größer oder gleich 3, so folgt vÎ Q. Gibt es jeweils nur eine Sequenz w und eine Sequenz x aus X mit v1 = w1 und v2 = x2, so müssen folgende Fälle unterschieden werden:

Sind w und x identisch, so folgt vÎ Q. Anderenfalls überprüfe man, ob es eine Merkmalsausprägung der ersten oder zweiten Koordinate gibt, die mehr als einmal im Datensatz vorkommt. Für den Fall, daß diese Ausprägung weder mit w2 noch mit x1 übereinstimmt, gilt vÎ Q. Sonst überprüfe man, ob es eine weitere Sequenz y in X gibt, die weder in der ersten noch in der zweiten Koordinate mit w oder x übereinstimmt. Ist dies der Fall, so folgt vÎ Q, anderenfalls endet der Algorithmus mit vÏ Q.

Der hier beschriebenen Algorithmus weist eine Komplexität von O((#L ) x (#X)) auf, wobei #X die Anzahl der unterschiedlichen Sequenzen in X sei. Für die Untersuchung einer Sequenz beliebiger Länge n errechnet sich dann in Kombination mit Lemma 1 ein Laufzeitverhalten des Algorithmus von O(n2 x (#L ) x (#X)).

 

Steinerbäume/ sparsamste (MP) Bäume in quasimedianen Netzwerken

 

Wie bereits erwähnt, ist es möglich, für drei Sequenzen mittels der quasimedianen Operation spezielle Lösungen des Steinerproblemes zu finden. Offen bleibt dabei die Frage, ob auch bei vorgegebenem beliebigem Alphabet L in einer durch die Sequenzen X Í L n erzeugten quasimedianen Algebra Q immer ein sparsamster (MP) Baum für eine gegebene Baumtopologie gefunden werden kann.

Für binäre Sequenzen, d. h. L = {0, 1}, ist es relativ einfach möglich, mit Hilfe des Algorithmus von Fitch (1971) zu einer vorgegebenen Baumtopologie mit gegebener Wurzel einen kürzesten Baum bzgl. der Anzahl der Merkmalsänderungen zu bestimmen. Hierzu werden alle gegebenen Sequenzen bis auf die Wurzelsequenz linear auf einer Geraden so angeordnet, daß sich je zwei vorgegebene Kanten nicht überschneiden (vgl. Fig. 2).

Dabei geht der Algorithmus so vor, daß Merkmalsweise die inneren Knoten des vorgelegten Baumes zunächst in einem ¢ Top-Down¢ -Verfahren, d. h. von den Endästen des Baumes zur Wurzel, mit den möglichen Ausprägungen belegt werden. Ist diese nicht eindeutig, wird eine Menge von Ausprägungen angenommen. Nicht eindeutig bestimmte innere Knoten werden dann in einem zweiten ¢ Bottom-Top¢ -Durchlauf durch die jeweils von unten angelieferte Ausprägung festgelegt.

Prinzipiell ist es nun möglich, dieses Verfahren, das ursprünglich nur für dichotome Bäume konzipiert war, auch für polytome Bäume und auch für andere Daten, z. B. DNS-Sequenzen, zu benutzen. Problematisch dabei sind aber zum einen die Bestimmung der möglichen Belegungen im Top-Down-Durchgang, als auch die eindeutige Festlegung der inneren Knoten im Bottom-Top-Durchgang, denn es kann vorkommen, daß ein innerer Knoten z. B. drei Ausprägungen {A, G, C} als Wahlmöglichkeiten für die Belegung zuläßt, das von unten angelieferte Merkmal aber nicht in dieser Menge enthalten ist (z. B. {T} ), so daß der innere Knoten mit drei Ausprägungen belegt werden könnte.

Während die Bestimmung der möglichen Belegungen des inneren Knotens noch relativ einfach durch die Majoritätsregel (¢ die am häufigsten von oben angelieferten Merkmale stehen zur Auswahl¢ ) gelöst werden kann, ist es für die eindeutige Festlegung der inneren Knoten notwendig, für jeden Knoten eine Prioritätsregel einzuführen, welche Merkmalsausprägung bei mehreren Alternativen zu wählen ist. Eine einfache Regelung diesbezüglich wäre: ¢ Die von links kommende Merkmalsausprägung wird in nicht eindeutigen Situationen übernommen¢ .

 

Beispiel:

Fig. 2: Künstliches Beispiel für den modifizierten Fitch-Algorithmus mit zusätzlichen Regeln wie oben angegeben. Untereinander stehende Merkmalsausprägungen geben alternative Wahlmöglichkeiten gemäß der Prioritätsregel an, wobei die Ausprägungen gemäß der Priorität von oben nach unten angeordnet sind. Durchgestrichene Merkmalsausprägungen wurden im Bottom-Top-Durchgang nicht gewählt. Die Prioritätsregel kommt in diesem Beispiel am zweiten Merkmal des inneren Knotens zum Tragen, von dem die Sequenz CG abzweigt.

 

Zu zeigen bleibt einerseits, daß der so modifizierte Fitch-Algorithmus weiterhin einen kürzesten Baum für die vorgegebene Baumtopologie liefert, zum anderen, daß alle inneren Knoten des so modifizierten Fitch-Algorithmus innerhalb der durch die Sequenzen X erzeugten quasimedianen Algebra enthalten sind.

Da die Einführung einer Prioritätsregel für jeden inneren Knoten nur dann zum Tragen kommt, wenn es mehrere äquivalente Wahlmöglichkeiten gibt, ist klar, daß diese Regel keinen Einfluß auf die Länge des Baumes hat. Zu zeigen bleibt, daß die Erweiterung des Fitch-Algorithmus auf polytome Bäume weiterhin eine sparsamste Realisierung liefert (vgl. Swofford et al. (1987)).

  

Die Belegung der inneren Knoten durch die Majoritätsregel läßt sich wie bei Fitch (1971) mathematisch ausdrücken durch die Intersunion Operation:

Für gegebene Mengen Æ ¹ L 1, ..., L k Í L , k ³ 2, ist durch

, wobei jÎ {1,...,k} die größte Zahl ist, für die die rechte Seite nicht die leere Menge ergibt, die Intersunion Operation definiert. Die Intersunion von L 1, ..., L n besteht also genau aus jenen Elementen, die in den n Mengen am häufigsten vertreten sind. Bei paarweise disjunkten Mengen L 1, ..., L k ergibt sich daraus j = 1 und der innere Knoten wird mit der Vereinigung der Mengen L 1, ..., L k belegt. Haben dagegen alle Mengen ein Element gemeinsam, so erhält man j = k und der innere Knoten wird mit dem Schnitt aller Mengen belegt.

Da somit ein innerer Knoten nur dann eine Merkmalsausprägung annehmen kann, wenn diese in den von oben angelieferten Sequenzen an dem entsprechenden Merkmal am häufigsten vorkommt, folgt, daß jede Wahl einer Merkmalsausprägung für die von oben angelieferten Sequenzen (Mengen) eine kürzeste Realisierung bzgl. der Anzahl der Zustandsänderungen der Merkmale liefert. Bei jeder Wahl einer Merkmalsausprägung, welche weniger häufig in den von oben angelieferten Sequenzen vorkommt, wäre mindestens eine Zustandsänderung mehr zur Verbindung der oben angelieferten Sequenzen nötig. Durch inverse Induktion bzgl. des Abstandes zur Wurzel folgt dies für alle inneren Knoten.

Im Bottom-Top Durchgang können alle inneren Knoten nur die im ersten Durchgang durch die Intersunion Operation festgelegten Ausprägungen eines Merkmals annehmen, d. h. für die von oben angelieferten Sequenzen ergibt die Wahl eines Elementes der durch die Intersunion Operation bestimmten Menge (kurz: Intersunion Menge) eine kürzeste Realisierung. Ist die Merkmalsausprägung der Wurzelsequenz ebenfalls in der Intersunion Menge enthalten, so wird der innere Knoten mit der eindeutig festgelegten optimalen Ausprägung belegt, und es ist nichts mehr zu zeigen. Ist dagegen die Ausprägung der Wurzelsequenz nicht in der Intersunion Menge enthalten, so wird gemäß der festgelegten Prioritätsreihenfolge ein Element dieser Menge gewählt. Zu zeigen bleibt, daß die Wahl der Ausprägung des Wurzelknotens keine kürzere Realisierung liefert. Dies folgt aufgrund der Eigenschaft der Intersunion Operation, für die von oben angelieferten Sequenzen stets eine kürzeste Realisierung zu garantieren. Bei der Wahl einer Ausprägung, die nicht in der Intersunion Menge enthalten ist, wächst dadurch die Anzahl der notwendigen Mutationen zur Verbindung der von oben angelieferten Sequenzen um mindestens eine Mutation an. Dieser Anstieg der Mutationszahl kann bestenfalls durch die Wahl der Ausprägung der Wurzelsequenz wieder ausgeglichen werden, und man hätte einen Baum mit der gleichen Länge gefunden, aber keinen kürzeren Baum. Durch Induktion bzgl. des Abstands zur Wurzelsequenz folgt dies für alle inneren Knoten, so daß der modifizierte Fitch Algorithmus für polytome Bäume in der Tat eine sparsamste Realisierung für jede vorgegebene Baumtopologie liefert. Q.e.d.

 

Lemma 3: Sei S eine nichtleere Teilmenge des Sequenzraumes L x L = L 2. Mit Si, iÎ {1, 2} sei die Projektionen von S auf das i-te Merkmal bezeichnet. Dann gilt: S ist genau dann eine Subalgebra von L 2, wenn S mit einer der folgenden Mengen übereinstimmt :

(1) {(l ,p (l )| l Î S1} für eine Bijektion p : S1 ® S2, oder

(2) {{l } x S1} È {S2 x {m }} = {x Î S½ x1 = l oder x2 = m } mit l Î S1 und m Î S2, oder

(3) S1x S2 mit #S1, #S2 ³ 2.

 

Beweis: ‘<=‘: Um zu zeigen, daß die Mengen des Typs (1), (2) und (3) Subalgebren von L 2 sind, muß gezeigt werden, daß die Mengen bzgl. der quasimedianen Operation abgeschlossen sind.

Dies folgt aufgrund der Eigenschaft der quasimedianen Operation. Für jede einzelne Koordinate hat diese die Eigenschaft eines dualen Diskriminators (vgl. Bandelt et al. (1994)), d. h. bei drei Wahlmöglichkeiten a, b und c wird genau dann die Ausprägung von b gewählt, wenn b und c gleich sind. In allen anderen Fällen wird die Ausprägung von a gewählt. Angewandt auf eine der obigen Mengen, folgt sofort, daß es nicht möglich ist, ein Element zu generieren, das nicht bereits in der entsprechenden Menge enthalten ist.

‘=>‘: Sei S eine Subalgebra von L 2 .

Für eine Sequenz l , m Î L sei

 

S1(m ) = {k |k Î L 1 mit (k ,m )Î S} und

S2(l ) = {n |n Î L 2 mit (l ,n )Î S}.

 

S hat genau dann die Form (1), wenn für alle l Î S1 und m Î S2 gilt:

#S2(l ) = #S1(m ) = 1.

Insofern nehme man an, daß es l , m Î L mit #S2(l ) + #S1(m ) ³ 3 gibt.

Gilt #S2(l ) ³ 2, so gibt es Elemente (l , b ) und (l , d ) aus S mit b ¹ d . Für jedes m Î S2 – {b , d } wähle man ein Element a Î L mit (a , m ) Î S. Dann ergibt sich

 

(l ,m ) = ((a , m ) (l , b ) (l , d ))Î S

 

und somit folgt { l } x S2 Í S. Analog ergibt sich S1 x { m } Í S im Falle #S1( m ) ³ 2.

Existieren zwei verschiedene Elemente k , l Î L , so daß S sowohl k x S2 als auch l x S2 enthält, so ergibt sich für jedes m Î S2 , daß gilt S1 x { m } Í S und somit S = S1 x S2. Ebenso folgt dies, wenn S mehr als eine Menge der Form S1 x { m } enthält. Somit hat S entweder die Form (2) oder (3). Q.e.d.

 

Satz 1: Der modifizierte Fitch-Algorithmus mit Prioritätsregel für die Festlegung der innerern Knoten einer vorgegebenen Baumtopologie der Sequenzen X generiert nur Knoten der durch X erzeugten quasimedianen Algebra Q.

 

Beweis: Aufgrund von Lemma 1 genügt es, dies für alle zwei-elementigen Teilmengen der Merkmale zu zeigen. Angenommen, es gäbe zwei Merkmale i und j aus der Menge der Merkmale {1,...,n}, für die obiger Satz nicht gelte.

Die durch die beiden Merkmale i und j generierte quasimediane Algebra S hat entweder die Form

 

(1) S=Li x Lj , oder

(2) S={(l ,p (l ))½ l Î Li, p (l )Î Lj} mit einer bijektive Abbildung p : Li ® Lj , oder

(3) S={{l } x Lj} È {Li x {m }}.

 

- Im Fall (1) werden die inneren Knoten mit Sequenzen der quasimedianen Algebra belegt, da diese alle Kombinationen der Ausprägungen der Ausgangsmerkmale enthält, und durch den Fitch-Algorithmus nur Merkmalsausprägungen gewählt werden können, die in den Ausgangssequenzen auftreten.

- Im Fall (2) gilt S=X. Da alle Ausgangssequenzen vom Typ (l , p (l )) sind, ist automatisch bei der Belegung der Sequenzen im Top-Down Durchgang die Ausprägung p (l ) als Wahlmöglichkeit zulässig, wenn auch l zulässig ist. Da die Abbildung p eine eindeutige Zuordnung der Elemente von Li und Lj darstellt, stimmt notwendigerweise auch die Prioritätsreihenfolge für die möglichen Belegungen der inneren Knoten überein, d. h. die bijektive Abbildung p : Li ® Lj stellt zusammen mit der vorgegebenen Prioritätsregel einen Ordnungs-isomorphismus dar. Im Top-Down Durchgang wird also jeder innere Knoten durch die Intersunion Operation mit einer Menge (M, p (M)) belegt, wobei M eine nichtleere Untermenge von Li ist.

Da auch der Wurzelknoten der vorgegebene Baumtopologie die Form (x, p (x)), xÎ Li hat, muß bei einer Übereinstimmung der ersten Merkmalsausprägung im Bottom-Top Durchgang auch immer eine Übereinstimmung in der zweiten Merkmalsausprägung vorliegen. Liegt keine Übereinstimmung mit der von unten angelieferten Merkmalsausprägung vor, so wird die Ausprägung gemäß der Priorität festgelegt. Da die Ausprägungen l und p (l ) aufgrund des Ordnungsisomorphismus immer die gleiche Priorität aufweisen, wird jeder innere Knoten automatisch ebenfalls mit einer Sequenz aus S belegt.

- Fall (3): Erster Teil des Beweises durch inverse Induktion nach dem Abstand zum Wurzelknoten. Jede Ausgangssequenz hat am ersten Merkmal die Ausprägung l oder am zweiten Merkmal die Ausprägung m . Da die möglichen Belegungen der inneren Knoten gemäß der Majoritätsregel bestimmt werden, weisen somit alle durch Ausgangssequenzen im Top-Down Durchgang durch die Intersunion Operation determinierten innere Knoten die Wahlmöglichkeit l am ersten Merkmal auf oder die Ausprägung m als mögliche Belegung der Ausprägung des zweiten Merkmals. Alle Belegungen haben also die Form

 

 

Angenommen, von den k ³ 2 oberen Nachbarn eines inneren Knotens u wurden i mit Mengen der Form a.), j mit der Form b.) und ki-j mit der Form c.) belegt. Die Intersunion dieser k Mengen liefert für die Belegung des Knotens u im Fall i > j eine Menge der Form a.), im Fall i<j eine Menge der Form b.) und eine Menge der Form c.), falls gilt i=j. Hieraus folgt die Induktion.

Zweiter Teil des Beweises durch Induktion nach dem Abstand zur Wurzel. Im ‚Bottom-Top‘-Durchgang wird zunächst der mit der Wurzel inzidierende Knoten festgelegt. Da die Wurzelsequenz ebenfalls entweder l als Ausprägung des ersten Merkmales aufweist und/oder m als Ausprägung des zweiten, liegt im Fall c.) der innere Knoten auch in S. Im Fall a.) wird stets unabhängig von der Wurzelsequenz l als Ausprägung des ersten Merkmals und im Fall b.) stets m als Ausprägung des zweiten gewählt, so daß auch in diesen Fällen der innere Knoten in S liegt. Durch Induktion bzgl. des Abstandes zur Wurzelsequenz folgt dies für alle innere Knoten

 

Im Widerspruch zur Annahme liegen also alle innere Knoten in S. Q.e.d

 

Da somit die Belegungen der durch den modifizierten Fitch-Algorithmus bestimmten inneren Knoten einer vorgegebenen Baumtopologie immer innerhalb der durch die Ausgangssequenzen X erzeugten quasimedianen Algebra liegen, folgt, daß innerhalb der quasimedianen Algebra auch für jede vorgegebene Baumtopologie mindestens eine sparsamste Lösung gefunden werden kann. Als Folgerung dessen ist demnach auch für jede Baumtopologie, die einen sparsamsten (MP) Baum der Ausgangssequenzen generiert, mindestens eine Lösung in der quasimedianen Algebra enthalten.

 

Netzwerk- Algorithmen

 

Median Joining Algorithmus (MJ)

 

Ein Verfahren zur heuristischen Näherungslösung des Steinerproblems ist der Median Joining Algorithmus (MJ) von Bandelt et al. (1999). Zur Konstruktion des Netzwerkes wird ein minimales aufspannendes Netzwerk als Startkonfiguration benutzt, definiert als Vereinigung aller kürzesten aufspannenden Bäume, welches relativ schnell und einfach durch eine naheliegende Modifikation des Algorithmus von Kruskal (1956) konstruiert werden kann. Innerhalb dieses Netzwerkes werden durch genau zwei Kanten verbundene Tripel von Sequenzen als Kandidaten für die Generierung von Medianvektoren ausgewählt und durch ein Zusatzkriterium weiter selektiert. Die Medianvektoren aus den gewählten Tripeln werden hinzugefügt, und der Algorithmus solange wiederholt, bis keine weiteren Medianvektoren mehr aufgrund der Auswahlkriterien generiert werden können.

Um lokale Optima umgehen zu können, ist ein Parameter e vorgesehen, der sowohl die Berechnung der Startkonfiguration (das minimale aufspannende Netzwerk) als auch das Auswahlkriterium der Medianvektoren beeinflußt, so daß bei hinreichend großer Einstellung des Parameters sogar das gesamte quasimediane Netzwerk generiert werden kann. Abgesehen von der deutlich anwachsenden Laufzeit beinhaltet dieses zwar aufgrund von Satz 1 für alle Baumtopologien eine sparsamste (MP) Realisierung, ist normalerweise aber zu groß und komplex für die Präsentation (bereits aus drei Sequenzen der Länge zwei können sechs Medianvektoren generiert werden).

Die wesentlichen Vorteile des Algorithmus bestehen sowohl in der Datenmenge, die verarbeitet werden kann (die implementierte Version 2.0 kann 14.000 verschiedene Sequenzen der maximalen Länge 500 verarbeiten), als auch in der Laufzeit: Für einen Datensatz mit 500 Sequenzen der Kontrollregion der mtDNS benötigt der Algorithmus üblicherweise bei Wahl e = 0 auf einem mit 120 Mhz getakteten PC mit Pentium Prozessor weniger als 1 Stunde. Durch eine wahlweise Gewichtung der Merkmale können zusätzlich Unterschiede der Mutationsraten einzelner Merkmale berücksichtigt werden.

 

Reduced Median Algorithmus (RM)

 

Der Reduced Median Algorithmus (RM) von Bandelt et al. (1995) wurde ebenfalls wie MJ für die Analyse intraspezifischer Daten in Form von binären Merkmalen konzipiert.

Der Algorithmus betrachtet dabei Projektionen des Datensatzes auf kleine Merkmalsmengen mit bestimmten Inkompatibilitäten und postuliert Parallelismen anhand der Frequenzen und der vorgegebenen Gewichtung. Dabei wird der Einfluß der Gewichtung mit einem Parameter r ³ 1 bestimmt.

Auch dieser Algorithmus geht sukzessive vor, d. h. in jedem Durchgang wird zur Auflösung von Inkompatibilitäten nur jeweils ein Merkmal in zwei Teilmerkmale aufgespalten. Zur Festlegung des aufzuspaltenden Merkmals geht der Algorithmus so vor, daß bei Vorliegen einer Inkompatibilität der Steinerbaum dieser Merkmale bzgl. der gegebenen Gewichtung ausgewählt wird, der zusätzlich durch die Frequenzen unterstützt wird. Aus der Menge der so zur Auswahl stehenden möglichen Aufsplittungen von Merkmalen zur Generierung von (lokalen) Steinerbäumen wird letztlich der ‚beste Merkmalssplit‘ aufgrund der Gewichtsverhältnisse und der Frequenzen selektiert und aufgesplittet. Das Verfahren wird solange wiederholt, bis aufgrund der Kriterien keine Aufsplittung von Merkmalen mehr möglich ist.

 

Vergleich von MJ mit RM

 

 

D(u, v) = d(u, v) - l (f(u) + f(v)),

 

wobei f(u) und f(v) die Frequenzen der Sequenzen u und v des Datensatzes X bezeichnen und l ein näher zu bestimmender hinreichend kleiner Skalierungsfaktor sei, so daß auch das abgewandelte Distanzmaß weiterhin nicht negativ ist.

 

Das Verfahren von Foulds et al. (1979) / greedy FHP

 

Das Verfahren von Foulds et al. (1979) wurde ebenso wie MJ und RM zur Berechnung einer heuristischen Lösung des Steinerbaum-Problems konzipiert. Dazu wird ebenso wie bei MJ zunächst mit Hilfe des Algorithmus von Kruskal (1956) das minimale aufspannende Netzwerk errechnet. Bei Foulds et al. wird aber der Algorithmus nur anhand von Beispielen beschrieben, so daß eine exakte eindeutig beschriebene Vorgehensweise nicht vorliegt. Einige wesentliche Grundzüge des Verfahrens wurden 1999 von Bandelt et al. zu Vergleichszwecken mit MJ zum greedy FHP-Verfahren zusammengefaßt:

Innerhalb des minimalen aufspannenden Netzwerkes werden ebenso wie bei MJ Tripel von Sequenzen (Originalsequenzen des Datensatz oder in einem vorhergehenden Durchlauf des Algorithmus generierte Medianvektoren) zur Generierung von Medianvektoren (‚Steiner points created by coalescing‘) ausgesucht, die in dem minimalen aufspannenden Netzwerk mit genau zwei Kanten verbunden werden können. Bei der Auswahl der Tripel, aus denen Medianvektoren generiert werden, greifen drei Kriterien:

 

 

Die aus den so ausgewählten Tripeln generierten Medianvektoren werden dem Datensatz hinzugefügt, und das Verfahren solange wiederholt, bis aufgrund der Auswahlkriterien keine Medianvektoren mehr generiert werden können (‘greedy FHP‘, vgl. Bandelt et al. (1999)). Sollte nach Abbruch des Verfahrens die zuvor errechnete untere Schranke für die Länge des Steinerbaumes noch nicht erreicht sein, so schlagen Foulds et al. weiterhin vor, in einem ‘Postprocessing‘-Verfahren, d. h. einem nachgeschalteten vom ersten Teil unabhängigen Verfahren, zu versuchen, die Länge des bisher errechneten Netzwerks zu verbessern. Dazu werden paarweise alle Sequenzen (auch zuvor generierte Medianvektoren) überprüft, ob die rechnerische Distanz im Netzwerk realisiert wurde. Ist dies nicht der Fall, so werden sukzessive weitere Medianvektoren hinzugefügt. Sollte sich bei Hinzunahme eines Medianvektors die Länge der kürzesten Bäume innerhalb des Netzwerks verkürzen, wird der Medianvektor beibehalten, anderenfalls wieder gelöscht. Dieses zusätzliche Verfahren ist aber im greedy FHP-Verfahren von Bandelt et al. 1999 nicht realisiert.

Vergleich von MJ mit dem greedy FHP-Verfahren

 

 

Vergleich von RM mit dem greedy FHP-Verfahren

 

 

Wahl eines Verfahrens

 

Aufgrund der praktischen Ergebnisse ist eine grundsätzliche Bevorzugung eines Algorithmus nicht gerechtfertigt, da die Güte der Ergebnisse abhängig vom jeweiligen Datensatz ist. In dem Tibetischen Datensatz von Torroni et al. (1994) liefert MJ mit e =0 alle MP-Bäume, aber auch drei Kanten, die in keinem MP-Baum vorkommen, während das greedy FHP-Netzwerk (FHP ohne Postprocessing) exakt aus der Vereinigung aller MP-Bäume besteht. In dem RM-Netzwerk mit Parameter r = 2 sind nicht alle kürzesten Bäume enthalten, so daß hier das greedy FHP Verfahren das beste Ergebniss erzielt (vgl. Fig.6 in Bandelt et al. (1999)). Bei Analyse des Afrikanischen Datensatzes von Chen et al. (1995) (vgl. Fig. 7) weist wiederum das RM-Netzwerk mit Parameter r = 2 eine Kante mehr auf als das MJ-Netzwerk mit e = 0 und enthält dadurch alle MP-Bäume. In dem MJ-Netzwerk mit e =1 ist diese Kante zwar auch enthalten, aber zusätzlich auch noch viele weitere Kanten, die in keinem MP-Baum vorkommen.

Während das RM-Verfahren bei Datensätzen mit jungen, genetisch gesehen eng miteinander verwandten Sequenzen in der Regel nicht so effektiv ist bzgl. der Anzahl der postulierten Parallelismen, haben sowohl MJ als auch FHP Probleme mit Datensätzen, die Sequenzen alter Populationen enthalten, da diese in der Regel lange Äste generieren, die in seltenen Fällen bereits aufgrund einer einzigen Mutation falsch zugeordnet werden.

Nach bisherigem Kenntnisstand sollten zumindest das RM-Netzwerk und das MJ-Netzwerk errechnet und verglichen werden, vorausgesetzt, die Daten erlauben eine derartige Vorgehensweise, d. h. es handelt sich im wesentlichen um binäre Merkmale.

Erfolgversprechend erscheint auch eine Kombination der Verfahren. Zunächst postuliert RM mit einer entsprechend hohen Wahl des Parameters r anhand der Gewichtung und der Frequenzen Parallelismen und nachfolgend errechnet MJ das Netzwerk des modifizierten Datensatzes (vgl. Rolf et al., vgl. Anhang B).

 

Star Contraction (SC) Algorithmus

 

Einerseits ist es wünschenswert, größere Datensätze zu analysieren, um genauere Ergebnisse zu erhalten bzw. bisherige Ergebnisse bestätigen zu können. Andererseits stellt sich dabei das Problem, einen großen Datensatz zu präsentieren. Eine mögliche Lösung dieses Problems besteht darin, bereits vor der eigentlichen Analyse durch einen Algorithmus (RM, MJ, FHP, Neighbour-Joining,...) sternförmige Sub-Phylogenien, die auf demographische Expansionen hinweisen, zu einer Sequenz zusammenzufassen. Einerseits wird dadurch die Datenmenge und somit auch die Laufzeit erheblich reduziert, andererseits wird bei einer hinreichend vorsichtigen Vorgehensweise der Informationsgehalt bzgl. der Phylogenie kaum beeinflußt.

Ziel dieses hierarchischen Vorgehens ist es, den Datensatz entsprechend so zu verkleinern, daß mit einem der Netzwerk-Algorithmen zunächst ein Übersichtsnetzwerk errechnet (und gezeichnet!) werden kann. Wenn dieses vorliegt, können einzelne Teile des Netzwerks näher analysiert werden. Der Sinn besteht also keinesfalls darin, Homoplasien innerhalb der Daten zu entfernen oder einen Datensatz um jeden Preis zu verkleinern. Ziel ist es nur, soviel Informationen wie möglich zu poolen, um die Analyse zu erleichtern.

Um dieses Grundkonzept adäquat umzusetzen, müssen insbesondere zwei Problematiken in Betracht gezogen werden:

 

  1. Es können Situationen auftreten, in denen Sequenztypen mehreren Sternen zugeordnet werden
  2. Der anzestrale Typ eines Sterns, d. h. das Sternzentrum, könnte nicht im Datensatz enthalten sein.

 

Um diesen Situationen gerecht zu werden, muß zunächst festgelegt werden, bis zu welcher Distanz D von einem potentiellen Sternzentrum überhaupt nach Sequenzen gesucht werden soll, die diesem Stern zugehören könnten, denn es macht wenig Sinn, alle Sequenzen eines Datensatzes als möglicherweise zugehörig zu betrachten. Als Distanzmaß d(x,y) sei hierbei der einfache ungewichtete evolutionäre Abstand, also die Anzahl der unterschiedlichen Merkmalsausprägungen der Sequenzen x = (x1,...,xn) und y = (y1,...yn) angenommen (Hamming-Distanz).

 

 

Mit n sei dabei die Sequenzlänge bezeichnet.

Aufgrund der Erfahrungen mit Daten des hypervariablen Segments I der mtDNA spricht vieles dafür, höchstens bis D =3 nach Sequenzen zu suchen, die einem potentiellen Stern zugeordnet werden. Zöge man nur zwei trennende Mutationen in Betracht, so werden Datensätze mit alten Expansionen, wie sie z. B. in Afrika vorkommen, unter Umständen nicht genug reduziert. Betrachtet man dagegen sogar Sequenzen, die sich durch vier oder mehr Mutationen unterscheiden, so würden in einem Datensatz mit einer jungen Population sehr viele Sequenzen zunächst falsch zugeordnet, wodurch nur der Rechen- und Speicheraufwand immens anstiege, und die Gefahr fehlerhafter Zuordnungen wüchse.

Wie bereits erwähnt, kann es aber sogar sinnvoll sein, D =5 anzunehmen, wenn ein Datensatz alte Expansionen enthält. Letztlich muß diese Vorgabe also vom Datensatz abhängig gemacht werden. Die Vorgabe D =3 kann nur als Anhaltspunkt dienen.

Um den Speicheraufwand gering zu halten, sei weiterhin als Mindestanforderung für einen Stern angenommen, daß wenigstens zwei Sequenzen im Datensatz existieren, die sich höchstens an einer Merkmalsausprägung vom potentiellen Zentrum unterscheiden (wie dies normalerweise im hypervariablen Segment I der mtDNS auch der Fall ist).

Der hier vorgestellte Star Contraction Algorithmus (SC) geht so vor, daß zunächst möglicherweise fehlende Sternzentren durch Medianvektoren rekonstruiert werden. Dabei wird die Suche auf die Tripel u,v,wÎ X von Sequenzen beschränkt, für die gilt

 

d(u,v) = 2 , d(u,w) £ D +1 und d(v,w) £ D +1.

 

Die erste Bedingung folgt aus der Annahme, daß bei Vorliegen eines Sterns mindestens zwei Sequenzen im Datensatz enthalten sind, die sich vom Sternzentrum an genau einer Merkmalsausprägung unterscheiden. Die anderen Bedingungen ergeben sich direkt aus der vorgegebenen Größe D des Suchfensters für sternförmige Strukturen.

In dem um die Medianvektoren erweiterten Datensatz werden dann die möglichen Sterne (Cliquen) berechnet. Dazu werden für jede Sequenz des erweiterten Datensatzes die Sequenzen des Datensatzes bestimmt, unterteilt in Distanzklassen 1,..., D , die sich an höchstens D Merkmalsausprägungen vom potentiellen Sternzentrum unterscheiden. Eine Clique mit Zentrum u setzt sich wie folgt zusammen:

 

 

In Hinblick auf die Speicherkapazität eines PC’s stellt sich nach der Berechnung der potentiellen Sterne die Frage, ab welcher Größe, d. h. ab welcher Anzahl der Einträge, eine Clique als potentieller Stern gelten soll. Alle Cliquen zu speichern, würde zu viel Speicherplatz erfordern und erscheint nicht sinnvoll. Ein einfaches (künstliches) Beispiel:

 

Fig. 3: Fiktives Netzwerk. Mutationen sind durch Querbalken angedeutet.

 

Obiger Baum enthält zwei Sterne mit Zentren bei Sequenz 5 und Sequenz 8. Sowohl die Clique mit Sequenz 3 als auch die mit Sequenz 11 als Zentrum erfüllen aber ebenfalls die Forderung, daß es mindestens zwei Sequenzen auf Abstand eins zum Zentrum gibt.

Der Unterschied zu den tatsächlichen Sternzentren liegt darin, daß die Verteilung der zugeordneten Sequenzen auf die Distanzklassen unterschiedlich ist.

Bei D =3 haben die Cliquen der Sequenzen 8 und 11 die Darstellungen:

 

 

C(8)

C(11)

Distanz 1

9; 11

8; 12

Distanz 2

5; 10; 12; 13

9

Distanz 3

3; 4; 7

5; 10; 13

 

Bei C(8) sind wesentlich mehr Sequenzen in der Distanzklasse zwei enthalten als in C(11), daher liegt es nahe, zur Selektion der Sterne die Anzahl der Sequenzen in höheren Distanzklassen ebenfalls mit einzubeziehen. Zu bedenken ist aber, daß die Anzahl der Fehlzuordnungen mit zunehmender Distanzklasse wächst, und es deshalb nicht sinnvoll ist, den Einträgen in höheren Klassen zuviel Beachtung zu schenken. Insofern bietet sich eine Kombination wie die folgende an:

 

mit z. B. t = 3,5.

Hierbei sei mit #Ck(u) die Kardinalität der Teilmenge Ck(u):={xÎ C(u)½ d(x,u)=k} der Clique von u bezeichnet, kÎ {1,...,D }. Da als Grundvoraussetzung bereits zwei Sequenzen in C1(u) vorhanden sein müssen, bedeutet dies, daß beispielsweise C2(u) mindestens drei Sequenzen enthalten muß oder C3(u) sechs Sequenzen.

Auch unter dieser Bedingung werden häufig noch Cliquen gebildet, die nicht als Sternzentren in Frage kommen, weil bisher nicht beachtet wurde, ob die Sequenzen einer ‚konkurrierenden‘ Clique besser zugeordnet werden können. Die Anzahl der in Betracht kommenden Cliquen wird aber bereits deutlich reduziert.

Deswegen müssen die Cliquen überprüft werden, ob Sequenztypen mehrfach zugeordnet wurden. Sollte dies der Fall sein, so bleibt der Eintrag nur in denjenigen Cliquen erhalten, in denen der Eintrag den geringsten Abstand zum Zentrum aufweist. Mehrfach zugeordnete Sequenzen dürfen allerdings dann (noch) nicht gelöscht werden, wenn das Cliquenzentrum eine Originalsequenz des Datensatzes ist, und die mehrfach zugeordnete Sequenz einen geringeren Abstand zu einem Cliquenzentrum aufweist, das aus einem zuvor generierten Medianvektor besteht. Durch das Löschen von Mehrfachzuordnungen können Medianvektoren als Cliquenzentren überflüssig werden, wenn nach der Beseitigung von Mehrfachzuordnungen nur noch ein Eintrag in der Clique des jeweiligen Medianvektors vorhanden ist. Das Löschen des Eintrags in der konkurrierenden Clique würde dazu führen, daß der gelöschte Sequenztyp gar keiner Clique zugeordnet wird (vgl. Fig. 4). Der dort schwarz gezeichnete Medianvektor erfüllt alle Bedingungen einer Clique und wird bei D =3 aufgrund von Sequenz 4, die Abstand vier zu den Sequenztypen 1 und 2 aufweist, generiert. Alle Sequenzen, die zunächst auch der Clique mit dem Medianvektor als Zentrum zugeordnet sind, werden mit Ausnahme des Sequenztyps 1 aufgrund der geringeren Distanzen zu den Cliquenzentren bei Sequenz 2 und 3 wieder entfernt, so daß der Medianvektor überflüssig wird. Wäre nun auch der Sequenztyp 1 aus der Clique bei 2 entfernt worden, so könnte dieser Sequenztyp erst im nachfolgenden Durchgang zugeordnet werden, obwohl er recht eindeutig zu dem Stern mit dem Zentrum 2 gehört.

 

Fig. 4: Fiktives Netzwerk. Mutationen sind durch Querbalken angedeutet, Medianvektoren sind durch schwarz ausgefüllte Kreise gekennzeichnet. Die Frequenzen sind proportional zu den Kreisgrößen.

 

In einem zweiten Durchgang werden somit zunächst die überflüssigen Medianvektoren gelöscht und der Durchgang ohne die Einschränkung bzgl. der Medianvektoren wiederholt.

Auch nach der bisherigen Vorgehensweise ist es möglich, daß es weiterhin Mehrfachzuordnungen gibt, nämlich genau dann, wenn ein Sequenztyp die gleiche Distanz zu mehreren Cliquenzentren hat. Diese Sequenzen können nicht eindeutig zugeordnet werden und müssen folglich aus allen Cliquen entfernt werden.

Da nunmehr die Cliquen eindeutig bestimmt sind, können alle in einer Clique befindlichen Sequenzen gelöscht und die Frequenzen auf das Cliquenzentrum aufaddiert werden.

Um die Kontraktion von Verzweigungspunkten des Netzwerkes zu verhindern, ist eine Nachuntersuchung nötig, die nach erfolgter Reduktion überprüft, ob Sequenzen zu Sternen zugeordnet wurden, die auf Kanten des resultierenden Netzwerkes liegen, bzw. davon abzweigen. Zu diesem Zweck muß zunächst das Netzwerk berechnet werden, bzw. die Kanten eines jeden Sternzentrums im resultierenden Netzwerk. Da die Laufzeit des SC-Algorithmus gering gehalten werden soll, wird allerdings darauf verzichtet, das gesamte Netzwerk mitsamt der Medianvektoren beispielsweise durch den MJ- Algorithmus zu berechnen. Statt dessen wird nur das minimale, eventuell auch das ‘e - relaxierte‘ minimale aufspannende Netzwerk errechnet, das ebenfalls im MJ Algorithmus als Startkonfiguration genutzt wird. Dieses erhält man durch Berechnung der e -zulässigen Kanten E eines jeden Sternzentrums.

Jede zu einem Stern mit Sternzentrum u zugeordnete Sequenz x muß jetzt dahingehend untersucht werden, ob es eine Kante e=uvÎ E gibt, auf der sich eine Merkmalsausprägung ändert, in der sich auch die beiden Sequenzen u und x unterscheiden. Ist dies der Fall, so wird die entsprechende Sequenz aus dem postulierten Stern entfernt, und die Frequenz des Sternzentrums wieder entsprechend verringert.

 

Algorithmus

 

Schritt 1: Zusammenfassung der Sequenzen auf Sequenztypen.

Schritt 2: Berechnung der Menge M aller Medianvektoren aus Tripeln u,v,wÎ X von Sequenzen des Datensatzes, für die gilt:

 

d(u,v) = 2 , d(u,w) £ D +1 und d(v,w) £ D +1

 

Alle errechneten Medianvektoren erhalten als Frequenz den Wert Null.

Schritt 3: Berechnung der Clique C(u) mit u als Cliquenzentrum für alle Sequenzen aus M È X. Mit #Ck(u) sei die Kardinalität der

Teilmenge Ck(u):={xÎ C(u)½ d(x,u)=k} der Clique C(u) bezeichnet, kÎ {1,...,D }. Gilt

 

 

mit etwa t = 3.5, so wird die Clique als potentieller Stern gespeichert, anderenfalls wird die Clique gelöscht.

Schritt 4: Löschen der Mehrfachzuordnungen von Sequenztypen in Cliquen, falls diese einen geringeren Abstand zu einem anderen

Cliquenzentrum aufweisen, und das Cliquenzentrum der anderen Clique entweder kein in Schritt 2 generierter Medianvektor ist, oder beide Cliquenzentren Medianvektoren aus M sind.

Schritt 5: Löschen der Cliquen mit Medianvektoren als Zentrum, die nur noch einen Eintrag aufweisen, da der Medianvektor nicht mehr

von den zur Clique zugeordneten Sequenzen unterstützt wird.

Schritt 6: Löschen der Mehrfachzuordnungen, falls die entsprechenden Sequenztypen einen geringeren Abstand zu einem anderen

Cliquenzentrum aufweisen.

Schritt 7: Löschen aller Einträge noch vorhandener mehrfach zugeordneter Sequenztypen, da eine eindeutige Zuweisung zu einer

Clique nicht möglich ist.

Schritt 8: Berechnung der e -zulässigen Kanten mit etwa e =0 der Sequenzen, die nicht zu Cliquen zugeordnet wurden. Überprüfung der

zugeordneten Sequenzen, ob diese auf den Kanten liegen oder davon abzweigen. Sollte dies der Fall sein, so werden die entsprechenden Sequenzen aus der entsprechenden Clique entnommen.

Schritt 9: Reduktion des Datensatzes um alle noch in den Cliquen vorhandenen Einträge. Diese werden gelöscht, und zu der Frequenz

des Cliquenzentrums wird die Summe der Frequenzen der zugeordneten Sequenztypen addiert.

Schritt 10: Erweiterung des Datensatzes um die Medianvektoren, die Zentrum einer nicht gelöschten Clique sind.

Schritt 11: Weiter bei Schritt 2, bis Abbruch durch Benutzer erfolgt oder keine weitere Reduktion mehr möglich ist.

 

Das gesamte Verfahren kann nach erfolgter Kontraktion natürlich noch beliebig oft auf den kontrahierten Datensatz angewandt werden, dabei ist aber zu beachten, daß jeder Durchgang des SC- Algorithmus zu einem Informationsverlust führen kann. Alternativwege im Netzwerk können nach Durchführung des Algorithmus fehlen.

In der Regel wird mit diesem Verfahren eine Reduzierung des Datensatzes um 40% - 50% der Datenmenge (Sequenzlänge multipliziert mit der Sequenzanzahl) erreicht, und offensichtlich fehlende Sternzentren werden bereits durch Medianvektoren ergänzt, so daß der nachfolgende Netzwerk-Algorithmus zur Errechnung der Phylogenie nur noch ca. 20% - 30% der Laufzeit für den gesamten Datensatz benötigt. Zusätzlich ist das Netzwerk eines derart zusammengefaßten Datensatzes sehr viel übersichtlicher als das entsprechende Netzwerk des gesamten Datensatzes. Eine detaillierte Analyse einzelner Teile des Netzwerkes kann nachträglich durchgeführt werden.

 

Vergleich mit Wills (1995)

 

Der Pruning- Algorithmus von Wills (1995) unterscheidet sich von diesem Ansatz ganz erheblich. Der Pruning- Prozess bei Wills bildet in einem ersten Schritt zunächst die Consensus- Sequenz des Datensatzes, d. h. es wird für jedes Merkmal der Majoritätszustand gebildet, und kontrahiert sukzessive den Datensatz in Richtung dieser Sequenz. Als ‘Pruning level’ wird dabei die Anzahl der trennenden Merkmale zwischen zwei Sequenzen bezeichnet. Es werden alle Paare von Sequenzen untersucht, die sich höchstens in einer dem Pruning level entsprechenden Anzahl Merkmale unterscheiden. Weicht eines der Merkmale von der Consensus- Sequenz ab, so wird es entfernt.

Während das Ziel von SC darin liegt, mit möglichst geringem Informationsverlust die Datenmenge auf das notwendigste Maß zu reduzieren, dient der Algorithmus von Wills zusätzlich zur Wurzelung eines Datensatzes. Problematisch bei der Vorgehensweise von Wills sind zum einen die Art, wie der Datensatz reduziert wird, zum anderen die Grundüberlegung bzgl. der Wurzelung.

Inhomogene Datensätze, bestehend aus unterschiedlich großen Stichproben verschiedener Populationen, werden bei einer derartigen Vorgehensweise die Bildung der Consensus-Sequenz erheblich beeinflussen, so daß die Wurzelung des Datensatzes anhand des von Wills vorgeschlagenen Verfahrens allenfalls in einigen Ausnahmefällen gute Ergebnisse erwarten läßt. Der Pruning-Algorithmus führt in einem inhomogenen Datensatz nur zu einer Kontraktion in Richtung der Consensus-Sequenz, also der überproportional vertretenen Population. Die gleiche Situation wird auch in einem Datensatz auftreten, in dem mehrere unterschiedlich starke Expansionen vorhanden sind, wie beispielsweise in Afrika (vgl. Chen et al. (1995)). Schwach besetzte, alte Sequenztypen bleiben nahezu unberücksichtigt, wogegen junge Expansionen wie in L2 und L3 überproportional in die Consensus- Sequenz eingehen. Daß in einem auf diese Art reduzierten Datensatz auch die Festlegung einer Wurzel problematisch ist, dürfte klar sein. Die Variation innerhalb des Datensatzes wird von Level zu Level geringer und ist bei hinreichend großer Wahl des Levels schließlich gleich Null, ebenso wie der Informationsgehalt des so generierten Datensatzes, der dann nur noch aus der Consensus- Sequenz besteht und automatisch mit der Wurzel übereinstimmt, egal mit welchem Verfahren diese bestimmt wird.

Ein weiteres erhebliches Problem des von Wills vorgeschlagenen Verfahrens besteht darin, daß selbst offensichtliche Parallelismen nicht als solche erkannt werden und daher der Datensatz falsch kontrahiert wird. In dem Datensatz von Chen et al. (1995) wird z. B. bei der Kontraktion übersehen, daß die RFLP-Schnittstelle 185l bei Sequenztyp 70 fast sicher einer Rückmutation unterliegt (der Gewinn der Schnittstelle relativ zu Anderson et al. (1981) wird durch eine Transversion charakterisiert und der nachfolgende Verlust durch eine Transition), so daß der Sequenztyp 70 durch den von Wills vorgeschlagenen Algorithmus nicht auf die Expansion bei Sequenz 71 (vgl. Fig. 5), sondern die gesamte Expansion bereits beim Pruning-Level 2 auf den Verzweigungsknoten zwischen den Sequenztypen 70 und 71 kontrahiert wird.

 

Beispiel: Chen et al. (1995)

 

Zur Demonstration der Verfahren wurde der Datensatz von Chen et al. (1995) ausgewählt. Er enthält RFLP-Daten, und zwar 79 Sequenzen der Länge 127, darunter 77 Sequenztypen bei Vernachlässigung der unzuverlässigen Schnittstelle 16517 HaeIII.

Wertet man den kompletten Datensatz mit dem MJ- Algorithmus mit dem Parameter e = 0 aus, so erhält man folgendes Netzwerk:

 

Fig. 5: MJ-Netzwerk mit e =0 des Datensatzes von Chen et al. (1995). Die Kreisflächen sind proportional zu den Frequenzen. Sequenzen, die die unzuverlässige Schnittstelle 16517 HaeIII haben, sind durch ein + gekennzeichnet. Die Pfeile geben an, bei welcher Sequenz die Schnittstelle des entsprechenden Restriktionsenzyms vorliegt. Mehrfach mutierte Merkmale sind unterstrichen. Der durch das schwarz ausgefüllte Dreieck gekennzeichnete Knoten repräsentiert die Consensus-Sequenz des Datensatzes. Der schwarz ausgefüllte runde Knoten entspricht der durch Midpoint rooting bestimmten Wurzel, der mit einem * ausgefüllte Knoten stellt die durch Primatensequenzen bestimmte Wurzelsequenz des Datensatzes dar. Der grau gepunktete Knoten repräsentiert die durch Midpoint rooting bestimmte Wurzel des durch den Pruning Algorithmus von Wills (1995) modifizierten Datensatzes mit Pruning Level 7. Die gestrichelte Linie stellt den einzigen Unterschied des MJ-Netzwerks zum RM-Netzwerk mit r=2 dar.

 

Der Pruning Algorithmus von Wills (1995) würde die alten langen Äste aus der Gruppe L1 bei entsprechender Wahl des Pruning Levels zugunsten der zwischen L2 und L3 liegenden Consensus-Sequenz reduzieren. Je größer der Pruning Level gewählt wird, desto drastischer wird die Fehlplazierung der Wurzel. Bei Pruning Level 7 liegt die durch Midpoint rooting bestimmte Wurzel schließlich innerhalb L2 (vgl. grau gepunkteten Kreis im Netzwerk) und steht damit im Widerspruch zu der Wurzelung durch Primaten-sequenzen und ‘Midpoint rooting’ (Röhl et al. in Vorbereitung, Chen et al. (1995)), sowie der Wurzelung von Sequenzen des Hypervariablen Segments I der mtDNS durch die Neanderthaler Sequenz durch Krings et al. (1997).

Wendet man den SC- Algorithmus mit dem Parameter D = 3 auf diesen Datensatz an, so erhält man nach einem Durchgang folgende Ausgabe:

Fig. 6: MJ- Netzwerk mit e =0 des durch den SC- Algorithmus nach einem Durchgang kontrahierten Datensatzes.

 

Der Datensatz wurde im ersten Durchgang des Algorithmus um 43 Ausgangssequenzen verringert und sowohl bei den Sequenzen 36 und 37 als auch bei 19, 20, 22 und 23 wurde jeweils ein Medianvektor generiert. Im reduzierten Datensatz sind also noch 38 Sequenzen mit 78 variablen RFLP- Schnittstellen enthalten.

Einigermaßen verwunderlich erscheint auf den ersten Blick, daß die Sequenzen 40 und 41 nicht der Expansion bei Sequenz 38 zugeordnet wurden. Der Grund hierfür liegt darin, daß Sequenz 41 sowohl zur Clique bei Sequenz 27 als auch zu Sequenz 38 jeweils Abstand drei aufweist. Das Merkmal 13065c ist wahrscheinlich zweimal mutiert, so daß die korrekte Distanz zu Sequenz 27 verfälscht wird.

Durch das Zusammenfügen der Sequenzen 4, 5, 6, 8, 9, 10 und 11 wird einer der beiden aneinanderhängenden Kreise aufgelöst. In diesem Fall ist dies durchaus sinnvoll, denn der Evolutionsweg innerhalb dieses Kreises wird bereits durch die Frequenzen so nahegelegt, daß auch der Cycle Pruning Algorithmus eine Auflösung in dieser Weise vorschlagen würde. Das resultierende Diagramm ist bereits wesentlich übersichtlicher, beinhaltet aber weiterhin alle wesentlichen Informationen des Gesamtdatensatzes.

Wird der SC- Algorithmus insgesamt dreimal auf den bereits reduzierten Datensatz angewandt, so ergibt sich das folgende Netzwerk:

 

Fig. 7: MJ- Netzwerk mit e =0 des durch den SC-Algorithmus nach drei Durchgängen kontrahierten Datensatzes.

 

Im zweiten Durchgang des Algorithmus wurden weitere 5 Sequenzen entfernt, so daß der resultierende Datensatz noch 33 Sequenzen der Länge 74 enthält. Ein weiterer Durchgang bewirkt noch, daß der aus den Sequenzen 36 und 37 gebildete Medianvektor nun der Expansion bei Sequenz 27 zugeordnet wird. Insgesamt wurde der Datensatz um 76.7% der Datenmenge (variable Merkmale multipliziert mit der Anzahl der Sequenztypen) reduziert.

Weitere Durchgänge haben keine Auswirkungen mehr. Das resultierende Netzwerk weist weiterhin alle tiefen Evolutionswege des Originals auf. Obwohl auch in der Haplogruppe L1c, die angesichts ihrer Diversität in diesem Datensatz nur durch die Sequenztypen 62 bis 69 repräsentiert wird, Expansionen stattgefunden haben, wurde dieser Bereich durch den SC-Algorithmus nicht kontrahiert, denn zum einen sind die Expansionen in diesem Datensatz nicht stark repräsentiert, so daß keine echten ‚Sterne‘ auftreten, zum anderen macht dieser Bereich des Netzwerks auch aufgrund der großen Anzahl trennender Merkmale untereinander wenig Probleme sowohl bei der Analyse als auch bei der Zeichnung. Eine Kontraktion dieser Sequenzen würde insofern auch bei der vorliegenden Datenlage keinen Sinn machen.

Als Nebeneffekt des SC-Algorithmus tritt bei diesem Datensatz eine Auflösung der beiden mit Sequenz 7 verbundenen Kreise auf, wodurch Parallelismen des Merkmales 12406 h/o postuliert werden, die auch bereits in dem Tibetischen Datensatz von Torroni et al. 1994 hochvariabel erscheint (vgl. Netzwerk in Bandelt et al. 1999).

 

Effektivitäts- und Laufzeitanalyse

 

Analysiert wurden Datensätze des hypervariablen Segments I im Bereich np16090 bis np16365 der mtDNS mit zufällig ausgewählten Sequenzen einer Datenbank, die insgesamt 2055 Sequenzen aus allen Teilen der Welt enthält, darunter 1279 Sequenztypen. Diese Datensätze wurden durch den SC-Algorithmus auf einem PC mit Pentium 120 Prozessor sowohl in einem Kontraktions-Schritt als auch in drei aufeinanderfolgenden Schritten zusammengefaßt.

 

Anzahl Sequenzen

Anzahl

Sequenz

typen

Variable Positionen

Zeit

(s)

Generierte Median-Vektoren

Anzahl

Durch-

Gänge

Anzahl

Sequenz-typen nach Kontraktion

Variable Positionen nach Kontraktion

Daten-reduktion

In %

200

173

123

92

10

1

134

98

38,29

250

210

124

126

11

1

157

100

39,71

300

244

139

235

10

1

168

97

51,60

350

280

129

458

12

1

209

107

38,09

400

301

134

542

6

1

236

115

32,71

450

358

148

748

11

1

274

122

36,91

500

387

150

943

15

1

289

123

38,76

550

412

155

1205

12

1

312

122

40,39

600

447

158

1703

13

1

331

129

39,54

650

475

151

2013

11

1

366

129

34,17

700

516

157

2690

11

1

389

139

33,26

 

Ergebnisse des SC-Algorithmus nach einem Schritt. Die Prozentangaben bzgl. der Datenreduktion beziehen sich auf die Datenmenge, die sich als Produkt der Sequenztypen mit der Anzahl der variablen Positionen errechnet.

Anzahl Sequenzen

Anzahl Sequenz-

Typen

Variable

Positionen

Zeit

(s)

Generierte Median-vektoren

Anzahl

Durch-gänge

Anzahl Sequenz-typen nach Kontraktion

Variable Positionen nach Kontraktion

Daten-reduktion

in %

200

173

123

144

10

3

111

90

53,05

250

210

124

264

11

3

142

98

46,56

300

244

139

427

13

3

149

92

59,29

350

280

129

875

13

3

198

106

41,89

400

301

134

965

10

3

193

102

51,19

450

358

148

1391

18

3

230

113

50,95

500

387

150

1928

23

3

243

115

51,86

550

412

155

2402

16

3

272

110

53,15

600

447

158

3329

15

3

281

122

51,46

650

475

151

4065

20

3

309

115

50,46

700

516

157

5483

15

3

339

127

46,86

 

Ergebnisse des SC-Algorithmus nach drei Schritten.

 

Bei der Anwendung des SC-Algorithmus auf Kontrollregion- Daten der mtDNS wird die Datenmenge (variable Positionen multipliziert mit der Anzahl der Sequenztypen) bereits im ersten Schritt in der Regel um fast 40 % reduziert. Die folgenden zwei Kontraktions-Schritte sind dagegen nicht mehr so effektiv. Es erfolgt im Durchschnitt noch eine weitere Reduktion um 10% der Datenmenge. Die Laufzeit wächst dagegen um ca. 100% relativ zu der Laufzeit für einen Durchgang.

Angesichts dieser Umstände und des Risikos, mit wiederholter Durchführung des Algorithmus einen Informationsverlust zu erleiden, sollte ein Datensatz durch den SC- Algorithmus zunächst nur in einem Schritt reduziert werden.

Mit einem Ansatz gemäß t » a x nb für die Laufzeit t des Algorithmus, wobei n die Datenmenge (variable Positionen multipliziert mit der Anzahl der Sequenztypen) bezeichne, erhält man mit der Methode der kleinsten Quadrate nach Anwendung des Logarithmus anhand der obigen Tabellen die Werte

a = 1.6 x 10-9 und b = 2.48

für die Daten nach einem Kontraktions-Schritt und

a = 0.9 x 10-9 und b = 2.59

für die Daten nach drei Kontraktions-Schritten.

Je nachdem, wie oft der Algorithmus angewendet wird, weist der Algorithmus demnach eine durchschnittliche Komplexität zwischen O(#X 2.5) und O(#X 2.6) für Daten des HVS I der mtDNS im Bereich np16090 bis np16365 auf.

Dies ist zwar schlechter als beispielsweise die Komplexität des MJ-Algorithmus von Bandelt et al. (1999), aber trotz der hohen Komplexität ist der SC-Algorithmus noch erheblich schneller bei Datensätzen mit bis zu 1000 Sequenzen. Z. B. hat der MJ-Algorithmus für den oben angegebenen Datensatz mit 700 zufällig ausgewählten Sequenzen des HVS I 6120 Sekunden benötigt bei e =0, der SC-Algorithmus dagegen nur 2.690 Sekunden. Zusätzlich verringert sich normalerweise bei der Durchführung des SC-Algorithmus die Homoplasie im Datensatz, da bereits offensichtliche Parallelmutationen postuliert werden. Die momentan implementierte Version des RM-Algorithmus hat (aufgrund schlechter Programmierung) z. B. eine Komplexität von O(n3) und wird somit erheblich vom SC-Algorithmus unterboten. Hier bietet sich in Hinblick auf die zu erwartende Rechenzeit bei großen Datensätzen eine vorgenommene Datenreduktion durch den SC-Algorithmus an.

 

Hochdimensionale Netzwerke / Cycle Pruning (CP) Algorithmus

 

Wie bereits in der Einleitung angedeutet, können Situationen auftreten, verursacht durch einen hohen Anteil Homoplasien (Rück- und Parallelmutationen, Sequenzierfehler), in denen das resultierende Netzwerk einen komplizierten und hochdimensionalen Torso hat, so daß eine Analyse des Datensatzes aufgrund der Problematik der visuellen Darstellung des Netzwerkes zumindest erheblich erschwert, eventuell auch kaum mehr möglich ist. Da der zuvor beschriebene SC-Algorithmus nicht dazu dient, Homoplasien zu reduzieren, sondern dies höchstens als Nebeneffekt auftreten kann, wird auch das Netzwerk des resultierenden Datensatzes i. a. noch sehr komplex sein. Um in derartigen Situationen trotzdem eine weitergehende Datenanalyse zu ermöglichen, wurde der nachfolgend beschriebene Cycle Pruning Algorithmus (CP) entwickelt.

Ausgehend von einem Netzwerk werden auf Basis eines vorgegebenen Kriteriums sukzessive Kanten entfernt, bis der gewünschte Grad der Auflösung vorliegt. Ein vergleichbarer Ansatz von Excoffier et al. (1994) benutzt hierbei alternativ als Kriterien die Frequenzen der Typen und die Populationszugehörigkeit, um zwischen phylogenetischen Pfaden im vorgegebenen Netzwerk zu wählen. Der CP-Algorithmus hingegen betrachtet ausschließlich die Frequenzen als Kriterium. Die Einbindung eines Geographie-Kriteriums ist momentan nicht vorgesehen, da die zugrunde liegenden Netzwerk-Algorithmen (MJ, RM) diese Information nicht mit ausgeben. Denkbar wäre höchstens eine zusätzliche nachträgliche Einbindung dieser Information per Hand.

In der Grundannahme zur Beurteilung der Güte von Kanten stimmen beide Algorithmen überein: Normalerweise wird ein häufig auftretender Sequenztyp eher Nachfahren haben als ein Typ, der nur selten gefunden wird. Insofern sollten also Kreise in der Regel so aufgelöst werden, daß Sequenztypen eher von einem hochfrequenten Typen abzweigen als von einem mit niedriger Frequenz. Bei Excoffier et al. werden allerdings nicht auftretende Rückmutationen in Betracht gezogen, so daß Kanten erhalten bleiben können, bei denen einer der Endknoten aufgrund einer offensichtlichen Rückmutation mit einem hochfrequenten Typen inzidiert.

Ebenso könnte der anzestrale Typ einer sternförmigen Expansion evtl. nahezu ausgestorben sein oder er wurde nicht gefunden und erst durch den Netzwerk-Algorithmus rekonstruiert (hier ist allerdings anzumerken, daß Excoffier et al. (1994) nur mit minimalen aufspannenden Netzwerken arbeiten, so daß letztere Problematik der Rekonstruktion von Sequenztypen hier nicht auftreten kann). Ein solcher Typ hätte nur eine geringe bzw. gar keine Frequenz, und in der Regel würden Kreise zu Ungunsten dieses Typs aufgelöst.

Um solche Problemsituationen zu kompensieren, berechnet der CP-Algorithmus zunächst für alle Knoten des Torsos die Attraktivität, indem die Frequenzen der Sequenztypen in der Umgebung betrachtet werden, wobei nicht nur die Anzahl trennender Mutationen, sondern auch die Gewichtung der Merkmale, mit einbezogen werden. Damit der Einfluß eines hochfrequenten Sequenztyps, der eine große gewichtete Distanz zu einem anderen Sequenztyp aufweist, nicht zu groß ist bei der Berechnung der Attraktivität, wird nur ein exponentiell mit dem gewichteten Abstand abnehmender Anteil der Frequenz addiert. Die Wahl des Faktors d Î (0,1), mit dem die Frequenz abnehmen soll, hängt hierbei vom Benutzer ab. Je nachdem, ob die Endäste viel oder wenig Einfluß auf die Kanten des Torsos ausüben sollen, muß der Faktor entweder nahe bei 1 oder nahe bei 0 gewählt werden. Der Wert d =0.5 ist nur als Richtlinie zu verstehen. Als Exponent des Faktors wird der gewichtete Abstand der Sequenztypen gewählt, so daß hoch gewichtete, trennende Merkmale auch bei der Berechnung der Attraktivität zu einer Art ‘Barriere’ werden.

Die Attraktivität a(v) eines Sequenztyps v errechnet sich dann durch

 

.

Dabei ist d(u,v) der gewichtete Abstand der Sequenzen u und v innerhalb des Netzwerkes, d. h. die Länge eines kürzesten Weges zwischen u und v.

Bei unitärer Gewichtung würde so die Frequenz f(u) eines Sequenztyps u, der sich in 4 Merkmalen von einem Sequenztyp v des Torsos unterscheidet, nur zu d 4 zu der Frequenz f(v) des Sequenztyps im Torso hinzuaddiert werden., d. h. die Attraktivität a(v) ergäbe sich als

a(v) = f(v)+d 4 f(u)

Zur Demonstration ein Beispiel:

 

Fig. 8: Fiktives Netzwerk. Medianvektoren sind durch schwarze Punkte gekennzeichnet und mit dem Zusatz ‘m‘ versehen. Die Kreisflächen sind proportional zu den Frequenzen.

 

Für obige Sequenzen errechnen sich mit d =0.5 die folgenden Attraktivitäten:

 

 

Frequenz

Attraktivität

Sequenztyp 1

1

8.9375

Sequenztyp 2

4

4.625

Sequenztyp 3

2

5.96875

Sequenztyp 4

2

5.96875

Sequenztyp 5

1

5.46875

Sequenztyp 6

4

6.96875

Sequenztyp 7

3

6.46875

Sequenztyp 8

2

7

Sequenztyp 9

1

5

Sequenztyp 10

1

4.4375

Sequenztyp 11

1

3.625

Sequenztyp 12

1

2.75

Sequenztyp m1

0

5.5

Sequenztyp m2

0

3.25

 

Für das Sternzentrum bei 1 errechnet sich eine hohe Attraktivität, obwohl Typ 1 selbst keine hohe Frequenz hat. Die hohe Frequenz von Sequenztyp 2 hingegen hat durch den großen Abstand (3 Mutationen) zum Torso nur geringen Einfluß.

Aufgrund der bisherigen Vorgehensweise können jetzt die Knoten bzgl. ihrer Glaubwürdigkeit als Verzweigungspunkte des Netzwerkes beurteilt werden, nicht jedoch die Kanten. Um ausgehend von den errechneten Attraktivitäten zu einem Kriterium für die Kantenbewertung zu gelangen, könnten zunächst die Attraktivitäten der Endknoten addiert werden. Bei einem Netzwerk, in dem nur gleich lange Kanten vorkommen, ist dieses Vorgehen auch ausreichend, nicht jedoch in einem Netzwerk mit unterschiedlichen (gewichteten) Kantenlängen. Hier muß mit einer Modifikation der Kantenbewertungen dafür gesorgt werden, daß eher lange als kurze Kanten reduziert werden. Eine Möglichkeit hierzu wäre die Einbindung dieser zusätzlichen Information durch abgewandelte Kantenlängen des Netzwerks gemäß

wobei a(u) und a(v) die errechneten Attraktivitäten der Knoten u und v des Netzwerks bezeichnen und l >0 ein näher zu bestimmender hinreichend kleiner Skalierungsfaktor sei, so daß die modifizierten Längen weiterhin positiv definit sind. Anders ausgedrückt: Einer Kante, die zwei Knoten mit einer hohen Attraktivität verbindet, wird ein vom Faktor l abhängiger ‚Rabatt‘ gegenüber der eigentlichen, rein auf Merkmalsänderungen basierenden Länge eingeräumt, aber nicht so viel, daß der resultierende Wert keinen Sinn mehr macht.

Ausgehend von den nun vorliegenden Kantenlängen kann jetzt ähnlich wie bei MJ das h -relaxierte minimale aufspannende Netzwerk errechnet werden, allerdings mit der Restriktion auf die bereits zuvor durch einen Netzwerk-Algorithmus errechneten Kanten, d. h. es können nur die Kanten in dem Netzwerk auftreten, die in dem zugrunde liegenden Netzwerk vorkommen, denn Sinn des CP-Algorithmus ist es, ein errechnetes Netzwerk aufgrund der zusätzlichen Frequenzinformation auszudünnen und nicht, ein völlig neues Netzwerk zu errechnen. Alle Kanten, die nicht in dem minimalen aufspannenden Netzwerk vorhanden sind, werden entfernt, wobei in dem resultierenden Netzwerk insbesondere darauf geachtet werden muß, ob Medianvektoren obsolet geworden sind.

Durch die Wahl des Parameters h kann der Grad der Ausdünnung bestimmt werden, indem mit einem entsprechend hohen Wert zunächst nur die ‚schlechtesten‘ Kanten entfernt werden und danach sukzessive der Wert von h herabgesetzt wird, bis der gewünschte Grad der Auflösung erreicht ist.

 

CP- Algorithmus

 

Schritt 1: Festlegung der Parameter d und l .

Schritt 2: Berechnung der Attraktivitäten a(v) für alle Knoten des Torsos gemäß

 

wobei f(u) die Frequenz des Knotens uÎ V bezeichnet und d(u,v) die Länge eines kürzesten Weges zwischen den Knoten u und v ist.

Schritt 3: Berechnung der modifizierten Kantenlängen gemäß

 

für alle Kanten (u,v) des Torsos.

Schritt 4: Festlegung des Parameters h und Berechnung des h -relaxierten minimalen aufspannenden Netzwerks.

Schritt 5: Entfernung aller Kanten des Netzwerks, die nicht im h -relaxierten minimalen aufspannenden Netzwerk enthalten sind.

Schritt 6: Überprüfung der Medianvektoren dahingehend, ob diese obsolet geworden sind oder nur genau zwei Kanten verbinden, so daß

diese zu einer zusamengefügt werden können.

Schritt 7: Berechnung des Torsos des neuen Netzwerks.

Schritt 8: Wiederholung des Verfahrens ab Schritt 2 bis das Netzwerk nur noch einen aufspannenden Baum darstellt oder bis der Benutzer

die Iterierung abbricht.

 

Verfahren von Excoffier et al. (1994)

 

Das 1994 von Excoffier et al. vorgeschlagene Verfahren wählt minimale aufspannende Bäume des Datensatzes als Lösungen aus, die bzgl. vorgegebener Kriterien optimal sind.

 

Dabei geht der Algorithmus im einzelnen so vor:

    1. Der Wert SSD(T) für die ‚total sum of squared deviations‘, der sich als Summe aller Produkte der relativen Häufigkeiten unterschiedlicher Sequenzen einer Population mit der im Baum gemessenen Distanz der Sequenzen ergibt.
    2. Der Wert SSD(WP) für die ‚within population sum of squared deviations‘, der sich ergibt, indem anlog zu 1.) für jede im Datensatz enthaltene Population der Wert SSD(T) berechnet wird und schließlich die Summe dieser Werte über alle Populationen gebildet wird.

Vergleich des CP-Algorithmus mit dem Verfahren von Excoffier et al. (1994)

 

Beide Algorithmen gehen zunächst davon aus, daß Kanten zwischen hochfrequenten Knoten wahrscheinlicher sind als solche, die nur durch wenige Sequenzen repräsentiert werden. Die Umsetzung bei Excoffier et al. ist allerdings insofern nicht in jedem Fall geeignet, da das vorgeschlagene Verfahren SSD(T) letztlich immer dazu führt, daß alle von einem hochfrequenter Typ abzweigenden Kanten (vgl. Sequenztyp 1 in Fig. 9) erhalten bleiben, auch wenn sich einige Sequenztypen vielleicht lokal wahrscheinlicher zuordnen lassen könnten und nur aufgrund einer Rückmutation eine Ähnlichkeit zu dem hochfrequenten Typ aufweisen. Hier geht der CP-Algorithmus lokaler vor: Ein hochfrequenter Typ hat zwar immer noch einen gewissen Einfluß auf weiter entfernt liegende Typen, aber dieser Einfluß nimmt exponentiell mit der Entfernung ab, so daß Rückmutationen eine Chance haben, als solche erkannt zu werden.

Bereits sinnvoller erscheint das zweite Kriterium SSD(WP), bei dem sich der errechnete Wert aus mehreren Einzelwerten zusammensetzt, nicht jedoch in der hier vorgeschlagenen Weise, da als Unterteilung des Datensatzes die Populationszugehörigkeit vorgeschlagen wird. Das hier benutzte Kriterium greift insofern nur zum Teil, da vielfach in Populationen mehrere Haplogruppen wiederum in Cluster unterteilt enthalten sind. Hier müßte insofern eine feinere Einteilung vorgenommen werden als die vorgeschlagene. Beim CP-Algorithmus kann dieses Problem nicht auftreten, da kein derartiges Kriterium vorgesehen ist.

Ein anderer wesentlicher Unterschied der Verfahren besteht aber auch darin, daß bei Excoffier et al. grundsätzlich eine Ausdünnung des Netzwerks bis auf einen minimalen aufspannenden Baum vorgesehen ist. Dies ist bei dem CP-Algorithmus nicht notwendig. Der Benutzer kann in gewissem Maße den Grad der Ausdünnung steuern und ist so in der Lage, problematische Kreise zu erhalten, deren Auflösung eher als Zufall anzusehen wäre. Das von Excoffier et al. vorgeschlagene Verfahren ist in dieser Art und Weise nur schwer modifizierbar, da die vorgeschlagenen Werte nur für Bäume berechnet werden können. Aus diesem Grund könnte man höchstens nachträglich wieder Bäume zu Netzwerken vereinigen.

Auch die Einbeziehung der Frequenzen in die Beurteilung der Kanten/ Bäume ist unterschiedlich. Der CP-Algorithmus errechnet zunächst Attraktivitäten der Sequenzen, d. h. eine Art ‚modifizierte‘ Frequenz, die zusätzlich berücksichtigt, wie die Frequenzverteilung um den entsprechenden Knoten herum aussieht. Durch dieses Vorgehen können zunächst alte Expansionsgründertypen rekonstruiert werden, d. h. Sternzentren mit einer sehr geringen Frequenz. Die Frequenzen werden bei dieser Berechnung allerdings nur exponentiell zum Abstand der Sequenzen abnehmend in die Berechnung einbezogen. Die Beurteilung der Kanten geht daher beim CP-Algorithmus lokaler vor als bei Excoffier et al., die die Verteilung der Frequenzen um die Endknoten einer Kante nicht berücksichtigen.

 

Beispiel: Excoffier et al. (1994)

 

Um die Ergebnisse des CP-Algorithmus mit dem von Excoffier et al. (1994) präsentierten Verfahren vergleichen zu können, wurde der dort analysierte Datensatz benutzt. Da der CP-Algorithmus ein bestehendes Netzwerk sukzessive ausdünnt, wurde als Startkonfiguration das MJ-Netzwerk mit e = 0 berechnet (vgl. Fig. 9).

Das von Excoffier et al. (1994) als Startkonfiguration benutzte und das MJ-Netzwerk unterscheiden sich ausschließlich in der Darstellungsart des Würfels bei Typ 11. Während MJ in diesem Bereich den vollständigen medianen Graphen generiert, werden bei Excoffier et al. drei ausgewählte Wege präsentiert.

 

Fig. 9: MJ Netzwerk mit Parameter e =0 des Datensatzes von Excoffier et al. (1994). Die Bezeichnung der Sequenztypen stimmt mit Excoffier et al. überein. Die Kreisflächen sind proportional zu den Frequenzen gewählt.

 

Aufgrund der wenigen benutzten Restriktionsenzyme und der dadurch auftretenden hohen Homoplasie bietet gerade dieser Datensatz eine gute Grundlage zum Testen des CP-Algorithmus. Da Excoffier et al. sowohl die Populationszugehörigkeit als auch die Frequenzen als Grundlagen zum Ausdünnen des Netzwerks benutzt haben, wird die Anwendung des CP-Algorithmus zusätzlich zeigen, inwieweit das von Excoffier et al. gezeigte Netzwerk von der impliziten Zugrundelegung insbesondere des zweiten Kriteriums abhängt.

Als Parameter wurden für den CP-Algorithmus in diesem Fall d =0.25, l =0.001, h 1=0.1, h 2=0.05, h 3=0.03, h 4=0.01, h 5=0.008, h 6=0.005, h 7=0.003, h 8=0.001 und h 9=0 gewählt, d. h. das Netzwerk wurde sukzessive in neun aufeinanderfolgenden Schritten ausgedünnt. Hierzu ist zu sagen, daß die Parameterwahl dem zugrundeliegenden Datensatz individuell angepaßt werden muß. In obigem Beispiel erscheint es aufgrund der geringen Auflösung sinnvoll, den Parameter d relativ klein zu wählen, so daß der Einfluß hoch frequenter Typen auf die Attraktivitäten weit entfernter Kanten sehr schnell abklingt. Der Faktor l muß in diesem Fall so gering gewählt werden, da das abgewandelte Distanzmaß sonst wegen der hohen Frequenz des Typs 1 nicht mehr positiv definit wäre. Die Werte für h wurden iterativ bestimmt, indem der Wert für h zunächst in Schrittweise herabgesetzt wurde, bis der Parameter wirkte. Aufgeführt sind hier nur die Parameterwerte, die eine Auswirkung auf das Netzwerk hatten. Bei einem anderen Netzwerk können insofern völlig andere Parameter sinnvoll sein. Bei Netzwerken mit langen Kanten und einer hohen Auflösung sollte der Parameter d eher höher angesetzt werden als bei einem Netzwerk wie das hier gezeigte mit einer geringen Auflösung. Die Wahl des Parameters l wiederum hängt von mehreren Kriterien ab:

 

 

Bei Wahl der oben angegebenen Werte der Parameter erhält man schließlich folgendes maximal ausgedünnte Netzwerk:

 

Fig. 10: Ergebnis des CP-Algorithmus mit Parametern d =0.25, l =0.001, h 1=0.1, h 2=0.05, h 3=0.03, h 4=0.01, h 5=0.008, h 6=0.005, h 7=0.003, h 8=0.001 und h 9=0 des Datensatzes von Excoffier et al. (1994). Die Bezeichnung der Sequenztypen stimmt mit Excoffier et al. überein. Die Kreisflächen sind proportional zu den Frequenzen gewählt.

 

Beim Vergleich des so erhaltenen Netzwerks mit den von Excoffier et al. präsentierten Netzwerken fällt auf, daß dieses Netzwerk ausschließlich Kanten enthält, die im SSD(T) oder im SSD(WP) optimierten minimalen aufspannenden Baum vorkommen, d. h. der CP-Algorithmus hat ein Netzwerk geliefert, das Aspekte beider Kriterien enthält. Entfernt man zusätzlich in dem übrig gebliebenen Kreis die Kante zwischen den Knoten 72 und 36, so erhält man einen minimalen aufspannnenden Baum und kann die entsprechenden Werte für SSD(T) und SSD(WP) wie bei Excoffier et al. angegeben berechnen. Man erhält SSD(WP)=321.937 und SSD(T)=411.213, also Werte, die kaum über den optimalen Werten von Excoffier et al. liegen. Somit hat der CP-Algorithmus für diesen Datensatz einen minimalen aufspannenden Baum geliefert, der bzgl. beider von Excoffier et al. vorgeschlagenen Kriterien nahezu optimal ist, verglichen mit den von Excoffier et al. präsentierten Maximalwerten für SSD(T)=2028.679 und SSD(WP)=1540.396. In der Addition der Werte (733.15) liegt der durch den CP-Algorithmus errechnete Baum im Bereich des bereits 1992 von Excoffier et al. vorgeschlagenen (732.07). Dabei liegt der Wert für SSD(T) des CP-Netzwerks um 0.607 unter und SSD(WP) um 1.687 über dem jeweiligen Wert des von Excoffier et al. vorgeschlagenen.

Auch ohne eine explizite Zugrundelegung der Populationszugehörigkeit kann also durchaus ein adäquates Ergebnis erzielt werden. Diese Zusatzinformation kann insofern bei einem nachträglichen Vergleich des Original-Netzwerkes mit dem durch den CP-Algorithmus errechneten, ausgedünnten Netzwerk eingehen, indem eventuelle Fehlzuordnungen bestimmt und per Hand korrigiert werden. In dem durch den CP-Algorithmus errechneten Netzwerk zum Beispiel wird die Verbindung des Typs 9 über eine Kante der Länge 1 mit Typ 6 statt über eine Kante der Länge 1 mit Typ 8 wie bei Excoffier et al. (1992) vorgeschlagen. Diese Kante bietet zwar eine Alternative, kann aber nicht unbedingt als bessere Lösung angesehen werden. Während Typ 6 zwar aus mehr Individuen als Typ 8 besteht, enthält Typ 8 im Gegensatz zu Typ 6 Tharu- Sequenzen aus dem Himalaya, die auch in Typ 9 enthalten sind. Hier wäre insofern eine alternative Plazierung zu erwägen.

 

Simulationsstudie

 

Bisherige Methoden zur Verifizierung der Datenmodelle bzgl. der phylogenetischen Analyse von Populationsdaten basieren im wesentlichen auf der Koaleszenztheorie (vgl. Wakeley et al. (1997)), in der, von einer Populationsstichprobe ausgehend, versucht wird, die Phylogenie rückwärts bis zum ersten gemeinsamen Vorfahren (‘most recent common ancestor’) zu verfolgen. Der triviale, aber aufwendige Weg, die Populationsevolution im Computer vorwärts zu simulieren, hat den Vorteil, Gründereffekte und Bottlenecks (mit dem einhergehenden Verlust von Sequenztypen) besser zu studieren.

 

Vorgehensweise

 

Zu diesem Zweck wird im folgenden eine Population mit einer Anfangsgröße von 500 Sequenzen (vier Sequenztypen mit je Frequenz 125) betrachtet, die verschiedenen Expansionen, Schrumpfungen und Stagnationen mit geringer Variabilität über einen Zeitraum von ca. 60.000 Jahren unterworfen wird, wobei als durchschnittliche Generationszeit ein Zeitraum von 25 Jahren angenommen wird. Bei der Simulation wird zugrunde gelegt, daß eine Frau maximal drei Töchter haben kann, die in das gebärfähige Alter kommen. Insgesamt steigt die Populationsgröße auf ungefähr 5000 Sequenzen an. Zusätzlich wird die von Forster et al. (1996) beschriebene Mutationsrate gewählt, die besagt, daß im hypervariablen Segment I der mtDNS im Bereich np16090 bis 16365 (Basenpaare nach Anderson et al. (1981) numeriert) im Mittel ca. alle 20180 Jahre eine Transition stattfindet.

Angenommen wird für die Simulation zunächst das ‚infinite sites‘- Modell, bei dem davon ausgegangen wird, daß Mutationsereignisse nie zweimal die gleiche Position treffen, d. h. die Anzahl der Nukleotid-Positionen ist unbeschränkt. Nach der Simulation wird jede Mutation per Zufallsgenerator mit einer von 275 Positionen identifiziert, wobei die Positionen in drei unterschiedliche Gruppen bzgl. der Wahrscheinlichkeit einer Transition eingeteilt werden. Die gewählten Wahrscheinlichkeiten und Gruppengrößen orientieren sich an der von Aris-Brosou et al. (1996) vorgeschlagenen Gamma-Verteilung mit Parameter a =0,4 für die Mutationswahrscheinlichkeit der Positionen. Als Grenzen für die maximalen Wahrscheinlichkeiten innerhalb der Gruppen wurden folgende Werte gewählt:

 

 

Über die so erhaltenen Gruppen wurden dann die Mittelwerte für die Wahrscheinlichkeiten gebildet. Da die von Aris Brosou et al. gewählte Gama-Verteilung sich nicht auf 275 Positionen des HVS I, sondern auf 300 Positionen bezieht, wurden die jeweiligen Bereiche mit dem Faktor 11/12 multipliziert und auf volle 5 Positionen gerundet.

 

Positionen 1-150

(konservativ)

Positionen 151-250 (durchschnittl.)

Positionen 251-275 (hochvariabel)

Wahrscheinlichkeit für eine Transition

1/1500

3/1000

3/125

 

Da bei dieser Indentifikation der Mutationen mehrere Merkmale mit derselben Position identifiziert werden können, tritt der Effekt von Rück- und Parallelmutationen auf. Andererseits ist aber die echte Evolution bekannt, so daß ein Vergleich der Ergebnisse der Algorithmen mit dem wahren Baum möglich ist

Aus dem resultierenden Datensatz wird schließlich ein ‚Sampling‘ von 200 Sequenzen per Zufallsgenerator gezogen und der phylogenetischen Analyse durch die zuvor beschriebenen Algorithmen RM, MJ, SC und CP unterzogen.

 

Für die Simulation über einen Zeitraum von 60.000 Jahren (2.400 Generationen) wurden diese Eckdaten vorgegeben:

 

 

Die Expansions- und Schrumpfungsphasen werden dabei gemäß einem Sinus-förmigen Verlauf im Intervall von [-p /2, p /2] simuliert. Zwischen den vorgegebenen Eckdaten ist bei der Simulation eine maximale Abweichung (zufallsabhängig) von 2% der Populationsgröße von Generation zu Generation erlaubt, d. h. es können zufallsbedingt auch während der Stagnationsphasen kleinere Expansionen und Schrumpfungen auftreten. So ist beispielsweise bei der hier beschriebenen Simulation eine leichte Expansion in den letzten 400 Generationen aufgetreten, so daß der Datensatz am Ende der Simulation 5371 Sequenzen enthielt. Im Laufe der Simulation sind außerdem zwei der ursprünglich vier Sequenztypen gänzlich ausgestorben, der erste bereits nach ca. 5.000 Jahren (200 Generationen) und der zweite nach ca. 20.000 Jahren (800 Generationen).

Aus dem resultierenden Datensatz wurden schließlich zweimal je 200 Sequenzen per Zufallsgenerator gezogen.

 

 

Datensatz 1:

 

 

Fig. 11: Originalbaum der zufällig gezogenen 200 Sequenzen des simulierten Datensatzes. Die Bezeichnung der Knoten entspricht der Nummer des ersten Repräsentanten dieser Sequenz innerhalb des Datensatzes. Die mutierten Positionen nach der Identifikation mit einer von 275 Positionen werden durch die Positionsnamen auf den Kanten repräsentiert. Die Kreisgrößen sind proportional zu den Frequenzen gewählt. Medianvektoren werden durch schwarz ausgefüllte Knoten repräsentiert. Die grau gefüllten Kreise entsprechen den Sequenztypen, die nach einem Durchlauf des SC-Algorithmus noch übrig sind. Pfeile weisen zu dem Sequenztyp, zu dem der entsprechende Typ durch den SC-Algorithmus zugeordnet wird. Die Sequenztypen 500 und 3373 entsprechen den Originaltypen, mit denen die Simulation gestartet wurde.

 

Anzumerken ist, daß die Erwartungshaltung falsch wäre, den obigen Baum exakt rekonstruieren zu können, den z. B. sind nach der Identifikation der Mutationen mit den 275 Positionen die Sequenztypen 1895 und 2388 als auch die Typen 2140 und 4815 identisch. Ebenso tritt auf dem Ast zwischen Typ 1190 und Typ 3954 eine Rückmutation an Position 264 auf, die aufgrund der Vorgehensweise, den Baum möglichst kurz, d. h. mit möglichst wenig Mutationen, zu konstruieren, kein Algorithmus in dieser Art errechnen, sondern vielmehr einen gemeinsamen Verzweigungspunkt (hier durch den Medianvektor dargestellt) postulieren wird.

Der SC-Algorithmus bewirkt bei diesem Datensatz eine Reduktion des Datenvolumens von 50 Sequenztypen mit 43 variablen Positionen auf 23 Sequenztypen mit 23 variablen Positionen, also eine Reduktion der Datenmenge um 75,4% (Sequenztypen multipliziert mit der Anzahl der variablen Positionen). Dabei wird kein Sequenztyp falsch zugeordnet und es werden bereits zwei Medianvektoren generiert. Im zweiten Durchgang wird der Datensatz um weitere 11 Sequenzen und 8 variable Positionen reduziert, der dritte Durchgang bewirkt noch eine weitere Reduktion um einen Sequenztyp und eine Position (nicht gezeigt).

Die Auswertung mit dem MJ-Algorithmus mit Parameter e =0 ergibt das folgende Netzwerk:

Fig. 12: MJ-Netzwerk des Datensatzes aus Fig. 11, errechnet mit Parameter e =0. Die Beschriftungen und die Darstellungsweise sind analog Fig. 11 gewählt.

 

Durch den MJ-Algorithmus werden einige Sequenztypen falsch verbunden und zwar Typ 353, der aufgrund einer Parallelmutation bei Position 260 fälschlich von Typ 3314 abzweigt und Typ 5110, der aufgrund der Parallelmutation bei Position 180 falsch von Typ 4043 abzweigt. Die letzte falsche Kante betrifft die Typen 3889 und 3050, die wegen der Mutation an Position 256 jetzt von Typ 5141 abzweigen statt direkt von Typ 3451.

Das vom MJ-Algorithmus errechnete Netzwerk enthält trotz der falschen Kanten einen Baum der exakt um drei Mutationen kürzer ist als der korrekte Baum, also die gleiche Länge hat wie der wahre Baum, abgesehen von den drei zuvor diskutierten Kante, die nicht rekonstruierbar sind. Außerdem zweigen die falsch zugeordneten Typen weiterhin von den richtigen Expansionszentren ab.

Im Vergleich dazu liefert der RM-Algorithmus mit Parameterwahl r=2 folgendes Ergebnis:

Fig. 13: RM-Netzwerk des Datensatzes aus Fig. 11, errechnet mit Parameter r=2. Die Beschriftungen und die Darstellungsweise sind analog Fig. 11 gewählt. Gestrichelte Kanten wurden durch den CP-Algorithmus mit Parameterwahl d =0.5, l =0.2, h 1=2.0, h 2=1.5, h 3=1.0, h 4=0.5, h 5=0.2, h 6=0 sukzessive entfernt. Fett gestrichelten Kanten geben die zusätzlichen Kanten der Kombination des RM Algorithmus mit Parameter r=2 und des MJ Algorithmus mit Parameter e =1 im Vergleich zum CP-Netzwerk an.

Das durch den RM-Algorithmus errechnete Netzwerk zeichnet sich durch einen vierdimensionalen Hyperwürfel aus, der das Netzwerk unübersichtlich macht und kaum glaubhafte Alternativwege aufweist. Im Gegensatz zu dem durch den MJ-Algorithmus errechneten und zum durch den CP-Algorithmus ausgedünnten Netzwerk beinhaltet allerdings das RM-Netzwerk den wahren Baum.

Der Parameter d wurde bei der Ausdünnung des RM-Netzwerks durch den CP-Algorithmus auf den Wert 0,5 gesetzt, da dieser Datensatz im Gegensatz zu dem Datensatz von Excoffier et al. (1994) eine hohe Auflösung aufweist. Alle anderen Parameter wurden der Wahl von d entsprechend angeglichen (vgl. S. 44)

Bei der Ausdünnung des RM-Netzwerks werden auch zwei korrekte Kanten entfernt, und zwar die Verbindung des Sequenztyps 2213 mit Typ 3373 über eine Kante der Länge zwei sowie die Verbindung des Typs 353 mit dem Medianvektor zwischen den Typen 3373 und 3065. Der Verlust dieser Kanten könnte allerdings umgangen werden, indem der CP-Algorithmus bereits früher abgebrochen wird, denn diese beiden Kanten werden erst in den letzten beiden Durchgängen entfernt. Mit Ausnahme der genannten Kanten weist das CP-Netzwerk ausschließlich Kanten des wahren Baumes auf.

Optimal arbeitet bei diesem Datensatz die Kombination des MJ- und des RM-Algorithmus (e =1 und r=2): Das resultierende Netzwerk ist einerseits deutlich weniger komplex als das RM-Netzwerk mit Parameter r=2, andererseits aber komplexer als das MJ-Netzwerk mit e =0 und beinhaltet so weiterhin den korrekten Baum, abgesehen von den drei Kanten, die nicht rekonstruiert werden können.

 

 

Datensatz 2:

 

Fig. 14: Originalbaum der zufällig gezogenen 200 Sequenzen des simulierten Datensatzes. Die Bezeichnung der Knoten entspricht der Nummer des ersten Repräsentanten dieser Sequenz innerhalb des Datensatzes. Die mutierten Positionen nach der Identifikation mit einer von 275 Positionen werden durch die Positionsnamen auf den Kanten repräsentiert. Die Kreisgrößen sind proportional zu den Frequenzen gewählt. Medianvektoren werden durch schwarz ausgefüllte Knoten repräsentiert. Die grau gefüllten Kreise entsprechen den Sequenztypen, die nach einem Durchlauf des SC-Algorithmus noch übrig sind. Pfeile weisen zu dem Sequenztyp, zu dem der entsprechende Typ durch den SC-Algorithmus zugeordnet wird. Die Sequenztypen 508 und 3573 entsprechen den Originaltypen, mit denen die Simulation gestartet wurde.

 

Auch in diesem Datensatz wird kein Algorithmus die Sequenztypen 1428 und 4310 sowie die Typen 2650 und 3907 als unterschiedliche Typen identifizieren können, da diese bei der Identifikation der Mutationen mit den 275 Positionen identisch geworden sind. Ebenso wird keiner der hier untersuchten Algorithmen die Rückmutation des Merkmales 264 als solche feststellen, sondern stets einen gemeinsamen Vorfahren postulieren (hier durch den in obigem Baum nicht notwendigen Medianvektor zwischen den Typen 3573 und 4303 gekennzeichnet).

Der SC-Algorithmus hat auch in diesem Datensatz ausschließlich korrekte Zuordnungen der Sequenzen vorgenommen und die Datenmenge bereits in einem Schritt von 42 Sequenztypen mit 39 variablen Positionen auf 19 Sequenztypen mit 22 variablen Positionen reduziert. Dies entspricht einer Reduktion der Datenmenge um 74.5%. In den beiden darauffolgenden Schritten wird der Datensatz noch um weitere 8 Sequenztypen und 6 Positionen reduziert (nicht gezeigt).

Die Auswertung mit dem MJ-Algorithmus mit Parameter e =0 ergibt das folgende Netzwerk:

 

Fig. 15: MJ-Netzwerk des Datensatzes aus Fig. 14, errechnet mit Parameter e =0. Die Beschriftungen und die Darstellungsweise sind analog Fig. 14 gewählt.

 

Ebenso wie bei Datensatz 1 werden auch hier durch den MJ-Algorithmus falsche Kanten postuliert und zwar exakt an den gleichen problematischen Stellen wie im Datensatz 1. Es fällt dabei auf, daß diese fehlerhaften Kanten ausschließlich die Sequenztypen betreffen, die eine größerer Distanz zum nächstgelegenen Sequenztyp aufweisen. Während alle Sequenztypen über Kanten des korrekten Baumes verbunden werden, die nur Abstand eins oder höchstens zwei zum nächstgelegenen Typen aufweisen, werden drei von fünf weiter entfernt gelegene Typen über falsche Kanten verbunden. Der Schluß liegt nahe, daß für Datensätze mit großen Astlängen der MJ-Algorithmus mit einer höheren Parameterwahl oder einer entsprechend anderen Gewichtung durchgeführt werden sollte, denn die Fehlerquote scheint bei großen Astlängen erheblich anzusteigen.

Das RM-Netzwerk mit Parameter r=2 sieht dagegen ebenso wie beim ersten Datensatz etwas komplexer aus:

 

 

Fig. 16: RM-Netzwerk des Datensatzes aus Fig. 14, errechnet mit Parameter r=2. Die Beschriftungen und die Darstellungsweise sind analog Fig. 14 gewählt. Dünn gestrichelte Kanten wurden durch den CP-Algorithmus mit Parameterwahl d =0.5, l =0.2, h 1=2.0, h 2=1.5, h 3=1.0, h 4=0.5, h 5=0.2, h 6=0 sukzessive entfernt.

 

Wie beim ersten Datensatz ist in dem RM-Netzwerk der wahre Baum enthalten, allerdings ist das Netzwerk auch wieder relativ komplex (wenn auch nicht in dem Maße wie beim ersten Datensatz), so daß auch hier die Anwendung des CP-Algorithmus sinnvoll erscheint. Die Ergebnisse des CP-Algorithmus sind mit denen des ersten Datensatzes vergleichbar. Auch hier werden zwei Typen falsch verbunden, wenn auch die korrekten Kanten erst bei einer sehr niedrigen Wahl des Parameters h entfernt werden.

 

Zusammenfassend kann gesagt werden, daß zum einen der MJ-Algorithmus im wesentlichen für die Netzwerkberechnung junger Datensätze benutzt, oder aber bei der phylogenetischen Analyse älterer Sequenzen mit einer hohen Auflösung ein höherer Wert für e gewählt werden sollte.

Der RM-Algorithmus dagegen hat bei dieser Simulation wesentlich komplexere Netzwerke errechnet, enthielt aber dafür auch in beiden Fällen den echten Baum. Die Kombination des RM-Algorithmus mit dem CP-Algorithmus erscheint aufgrund der Ergebnisse dieser Simuationsstudie erfolgversprechend. Allerdings sollte ein komplexes Netzwerk nicht wie hier maximal ausgedünnt werden, sondern man sollte noch einige unsichere Kreise übrig lassen und diese später per Hand untersuchen. In den obigen Beispielen würde allerdings wahrscheinlich auch per Hand der Fehler gemacht, keine Parallelmutation der Position 180 anzunehmen, da diese ja als einigermaßen konservative Position eingestuft wurde und es insofern naheliegt, die Kreise zu Gunsten der Position 180 aufzulösen.

Optimal zumindest in dieser Simulationsstudie ist die Hintereinanderausführung zunächst des SC-Algorithmus und danach des MJ- oder auch des RM-Algorithmus. Da durch den SC-Algorithmus bei der Kontraktion in obigen Datensätzen bereits viele Parallelmutationen postuliert werden, ist die Berechnung des Netzwerks anhand des kontrahierten Datensatzes unproblematisch. Bemerkenswert ist, daß der SC-Algorithmus trotz hoher Homoplasie der Daten die Sequenztypen absolut korrekt zuordnet und sogar eine Kontraktion der Datenmenge von mehr als 70% bereits im ersten Durchgang bewirkt.

Auch die Hintereinanderausführung des MJ- und des RM-Algorithmus mit den Parametern e =1 und r=2 hat in dieser Studie ein adäquates Ergebnis geliefert, so daß auch diese Kombination der Algorithmen in Erwägung gezogen werden sollte.

 

Anhang A: NETWORK 2.0

 

Das Programmpaket NETWORK 2.0 wurde komplett in Borland Pascal 7.0 geschrieben und ist für die Analyse intraspezifischer, rekombinationsfreier Daten, wie mtDNS (in Form von Punktmutationen, Aminosäure- Sequenzen und RFLP- Schnittstellen) oder Y-chromo-somalen STR‘s bzw. Punktmutationen konzipiert. Es kann auf jedem DOS-Rechner ab 80386er Prozessor installiert werden. Allerdings sollten mindestens 1 MB Hauptspeicher und ca. 3 MB auf der Festplatte verfügbar sein.

Durch den Aufruf ‘network’ innerhalb des entsprechenden Verzeichnisses auf der Festplatte wird eine Benutzeroberfläche gestartet, von der aus alle implementierten Werkzeuge und Algorithmen aufgerufen werden können. Das komplette Programmpaket ist mit Ausnahme einiger Optionen des Zeichenalgorithmus ausschließlich mit der Tastatur zu bedienen. Die wesentlichen Funktionstasten sind die Cursor-Tasten und die Escape- Taste.

In NETWORK sind Werkzeuge zur Generierung der Datensätze enthalten, die die benötigten Daten abfragen und in einem entsprechendem Format speichern. Insofern sollte also ein Benutzer nur dann einen Datensatz per Hand generieren, wenn er einen bestehenden Datensatz aus einem anderen (unbekannten) Format konvertieren möchte.

 

Datenformate

 

a.) MJ- Format, konzipiert für den MJ- Algorithmus, ab Version 2.0.

 

 

Es ist darauf zu achten, daß der Datensatz mit der letzten Gewichtszeile abschließt. Danach dürfen keine weiteren Zeilen (auch keine Leerzeilen!) angefügt werden, sonst treten Einlesefehler auf.

 

 

Beispiel : Datensatz von Saitou et al. (1997), Aminosäuresequenzen.

 

1

1

1

1

1

1

1

1

1

1

2

2

2

2

2

2

2

2

2

5

5

5

6

6

7

7

9

9

9

1

1

1

3

4

6

6

7

8

3

6

7

3

9

4

6

5

7

8

0

4

6

5

0

6

8

6

3

1

 

 

 

 

 

 

T

P

A

T

Q

E

R

F

E

R

V

M

F

G

S

L

G

E

A

 

 

1

2

 

 

 

 

 

 

T

L

A

T

Q

E

R

F

E

R

V

M

F

G

S

L

G

E

A

 

 

1

3

 

 

 

 

 

 

T

L

A

T

Q

E

R

F

E

R

V

M

F

G

S

L

G

E

A

 

 

1

4

 

 

 

 

 

 

T

P

A

T

Q

E

R

F

E

R

V

M

A

G

S

L

G

E

A

 

 

1

5

 

 

 

 

 

 

T

P

A

T

Q

E

R

F

E

R

V

M

F

G

S

L

G

E

A

 

 

1

6

 

 

 

 

 

 

T

L

A

T

Q

E

R

F

E

R

V

M

F

G

S

L

A

E

A

 

 

1

7

 

 

 

 

 

 

T

P

A

T

Q

E

G

F

E

R

V

M

F

S

S

M

A

E

A

 

 

1

8

 

 

 

 

 

 

T

P

A

T

Q

E

G

F

E

R

V

M

F

G

S

M

A

E

A

 

 

1

9

 

 

 

 

 

 

T

P

A

T

Q

E

G

F

E

R

V

M

F

S

S

M

A

E

A

 

 

1

1

0

 

 

 

 

 

T

P

A

T

Q

E

G

F

E

R

V

M

F

G

S

L

R

E

A

 

 

1

1

1

 

 

 

 

 

T

P

A

T

Q

E

R

F

Q

R

V

M

F

G

S

L

G

E

A

 

 

1

1

2

 

 

 

 

 

T

P

P

T

Q

E

R

F

Q

R

V

M

F

G

S

L

G

E

A

 

 

1

1

3

 

 

 

 

 

T

L

A

T

Q

E

R

F

Q

R

V

M

F

G

S

L

G

E

A

 

 

1

1

4

 

 

 

 

 

T

P

A

T

Q

E

R

F

Q

R

V

M

F

G

S

L

G

E

A

 

 

1

1

5

 

 

 

 

 

T

P

A

T

Q

E

R

F

E

R

V

M

F

G

S

M

A

E

A

 

 

1

1

6

 

 

 

 

 

T

P

A

T

Q

E

R

F

E

R

V

M

F

G

S

M

A

E

A

 

 

1

1

7

 

 

 

 

 

T

P

A

T

Q

E

R

F

E

R

V

M

F

G

S

M

A

L

A

 

 

1

1

8

 

 

 

 

 

T

P

A

T

Q

G

R

F

E

R

V

M

F

G

T

L

G

E

T

 

 

1

1

9

 

 

 

 

 

T

P

A

T

Q

G

R

F

E

R

V

M

F

G

T

L

G

E

A

 

 

1

2

0

 

 

 

 

 

S

P

A

T

L

G

L

L

E

R

V

I

F

G

T

L

G

E

A

 

 

1

2

1

 

 

 

 

 

T

P

A

A

Q

G

R

F

E

Q

V

M

F

A

S

L

G

E

A

 

 

1

2

2

 

 

 

 

 

T

P

A

A

Q

G

R

F

E

R

A

M

F

A

S

L

G

E

A

 

 

1

2

3

 

 

 

 

 

T

P

A

A

Q

G

R

F

E

R

A

M

F

A

S

M

A

E

A

 

 

1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0

1

0

1

0

1

0

1

0

1

0

1

0

1

0

1

0

1

0

1

0

1

0

1

.

.

.

.

.

 

b.) RFLP Format . Ursprünglich kreiert von A. Torroni, Università di Roma.

 

Jede Sequenz wird in drei Zeilen angegeben.

 

Bei diesem Datenformat muß beachtet werden, daß jede Sequenz in exakt drei Zeilen aufgeschrieben wird, d. h. die Gesamtzahl der Zeilen des Datensatzes muß ein Vielfaches von drei sein. Wichtig für das korrekte Einlesen der Daten ist außerdem die Beachtung der Trennung bei der Angabe mehrerer Schnittstellen durch ein Leerzeichen.

 

Beispiel:

 

9 (Name)

+16517e (Mutationen relativ zur Consensus-Sequenz)

1 (Frequenz)

13

-663e +16517e +9bpdel

1

33

-663e +10394c +10397a +13259o +13633e

1

97

-9052n -11403g

1

117

-663e -8572e +8920e +16517e +9bpdel

1

118

-663e +10394c +16517e +9bpdel

1

119

-663e -8221c +10394c +10397a +11321a +13259o +16517e

1

 

Tip : Um später problemlos mehrere Datensätze zu einem großen Datensatz vereinigen zu können, sollten die Schnittstellen in allen Datensätzen gleich angegeben werden, d. h. es ist nicht sinnvoll, in einem Datensatz den Gewinn einer Schnittstelle (relativ zur Sequenz von Anderson et al. (1981)) zu schlüsseln und in einem anderen Datensatz den Verlust derselben.

 

c.) Y-chromosomales Format, konzipiert für STR‘s.

 

Analog zum RFLP-Format wird jeder Sequenztyp auf drei Zeilen angegeben.

 

Ab der vierten Zeile werden die Sequenzen angeben:

 

 

Jede Sequenz muß genau drei Zeilen einnehmen. Punktmutationen können dadurch eingebettet werden, daß man diese als Repeat mit der Repeatanzahl 0 angibt und nur bei den entsprechenden Sequenzen, die eine Mutation aufweisen, den Wert 1 einträgt.

 

Beispiel:

 

DY390,DYS19,DY389,DY392,DY393,DY391,DY156.................(Loci- Bezeichnungen)

(Leerzeile)

(Leerzeile)

abo43 ................................................................................................(Name)

12,15,2,13,13,10,99 ........................................................................(Repeat- Angaben)

1..........................................................................................................(Frequenz)

abo51

12,15,2,13,13,10,99

1

ger1,6

12,16,4,11,13,10,7

1

ge1,54

10,14,2,13,13,9,7

1

ge1,57

11,14,2,13,13,11,7

1

 

MJ- Algorithmus (Bandelt et al. (1999))

 

Der MJ- Algorithmus ist bis auf wenige Ausnahmen genauso implementiert, wie in der Veröffentlichung beschrieben. Die maximale Sequenzlänge wurde von Länge 250 auf Länge 500 erweitert, so daß die in Bandelt et al. (1999) angegebenen Laufzeiten nicht mehr aktuell sind, wenngleich auch die durchschnittliche Komplexität von O(n2.2) für HVS I Sequenzen davon nicht beeinflußt wird. Die maximale Sequenzanzahl beträgt nun 14.000 verschiedene Sequenzen statt bisher 9.000, das 1979 von Foulds et al. vorgeschlagene Selektionskriterium kann alternativ zum publizierten Kriterium ‘connection cost’ gewählt werden, und der Satz der zulässigen Zeichen für DNS-Sequenzen wurde für nicht (oder nicht eindeutig) bestimmte Merkmalsausprägungen gemäß der in Anderson et al. (1981) veröffentlichten Abkürzungen erweitert. Zusätzlich kann optional eine Selektion der Sequenzen durchgeführt werden, die mindestens zweimal vorkommen bzw. eine höhere Frequenz als Eins aufweisen. Da die Wahrscheinlichkeit eher gering einzuschätzen ist, bei der Sequenzierung den gleichen Fehler mehrfach zu machen, kann man mit dieser Option versuchen, Homoplasien zu entfernen, die durch fehlerhafte Daten verursacht werden. Der Algorithmus arbeitet je nach Verfügbarkeit des Speichers entweder mit einer Distanzmatrix oder errechnet die jeweils benötigte Distanz neu. Da die zuletzt beschriebene Vorgehensweise ca. zehnmal langsamer ist als die Speicherung der Distanzen in einer Matrix, sollte möglichst viel Speicher zur Verfügung stehen. Der Algorithmus prüft nach jedem Schritt, ob noch genügend Speicher für die Distanzmatrix vorhanden ist. Bei Erreichen der Grenze tritt dann ein Sprung in der Laufzeit auf.

Bei der Benutzung des Algorithmus ist zu beachten, daß die Wirkung der e - Schranke direkt von der Gewichtung der Merkmale abhängt. Beträgt das kleinste Gewicht eines Merkmals z. B. 5, so hat der Wert e =4 keine andere Wirkung als e =0. Bei nahezu unitärer Gewichtung (fast alle Merkmale haben Gewicht 1) hingegen würde die Wahl e =4 wahrscheinlich bereits die Berechnung des kompletten quasimedianen Netzwerkes (oder zumindest eines nicht geringen Teiles davon) verursachen. Auswirkungen auf die Berechnung des Netzwerkes haben insofern nur Werte, die den Differenzen innerhalb der Distanzmatrix angepaßt werden.

Möchte man eines oder mehrere Merkmale ganz aus der Berechnung ausschließen, so kann man dies durch eine Null-Gewichtung erreichen. Merkmale, deren Gewichte Null betragen, gehen nicht in die Berechnung des Netzwerkes ein.

Die Ausgabe des Netzwerks kann wahlweise über einen Drucker, den Bildschirm oder in eine Datei erfolgen und enthält Listen der gleichen Sequenzen und der Kanten inklusive der auf der jeweiligen Kante auftretenden Merkmalsänderungen sowie eines Kennzeichens, ob die Kante innerhalb des Torsos liegt. Bei der Ausgabe in eine Datei werden zusätzlich die Frequenzen und die Gewichtung ausgegeben. Weiterhin enthält der Dateikopf die Einstellungen der Parameter.

Grundsätzlich ist eine Ausgabe in eine Datei zu empfehlen, da diese später immer noch ausgedruckt werden kann, und die Ausgabedatei des MJ-Algorithmus als Eingabedatei für das Graphikprogramm dient.

 

RM-Algorithmus (Bandelt et al. (1995))

 

Der Reduced Median (RM) Algorithmus wurde bereits 1996 von Duc Luu Minh an der Universität Hamburg in ein Programm umgesetzt. Aufgrund mehrerer Fehler wurde diese Version im Dezember 1998 neu programmiert und entspricht jetzt dem publizierten Algorithmus mit mehreren Verbesserungen:

 

 

Auch hier werden wie bei dem MJ- Algorithmus Merkmale vor der Berechnung aussortiert, deren Gewichte Null betragen. Die Ausgabedatei mit der Endung ‘out’ entspricht genau dem Format, welches sowohl das Graphikprogramm zur Zeichnung der Netzwerke, als auch der Cycle Pruning Algorithmus zur Ausdünnung eines Netzwerks benötigt. Die zweite Datei mit der Endung ‘mat’ enthält die veränderte Eingabematrix der binären Sequenzen in MJ-Format. Eine Kombination der beiden Algorithmen ist daher relativ einfach möglich. Der RM-Algorithmus kann sowohl das Y-chromosomale Format als auch das RFLP-Format problemlos verarbeiten. Das MJ-Format ist hingegen nur zulässig, wenn der Datensatz binäre Daten in Form von 0/1-Sequenzen enthält.

 

Graphikprogramm

 

Das Graphikprogramm benutzt genau das Ausgabeformat des MJ- Algorithmus als Eingabeformat. Um alle Optionen des Programmes nutzen zu können, sollte eine IBM- kompatible Maus installiert sein. Bisher haben z. B. Mäuse der Firmen Vobis, Microsoft und Logitech keinerlei Probleme gemacht. Voraussetzung ist aber die korrekte Installation der Maus. Während MS-Windows zur Not über eigene Maustreiber verfügt, ist dies bei NETWORK nicht der Fall. Sicherzustellen ist, daß die mit der Maus gelieferten Maustreiber ‘mouse.sys’ in der ‘config.sys’ und ‘mouse.com’ in der ‘autoexec.bat’ aufgerufen werden (die hier angegebenen Namen können eventuell bei anderen Herstellern abweichen). Ist dies nicht gewährleistet, werden Probleme bei der Benutzung auftreten.

 

Das Programm läuft gemäß folgendem Schema:

 

 

Als Optionen für die Darstellung eines fertigen Diagramms sind Vergrößerung, Verkleinerung, Verschiebung des Bildmittelpunktes, Ein- und Ausblenden der Sequenznamen und Medianvektoren sowie die Beschränkung des Netzwerkes auf den Torso vorgesehen. Der Ausdruck der Netzwerke kann wahlweise direkt auf den Drucker oder über MS-Windows erfolgen, wobei der direkte Ausdruck nur auf Nadeldruckern und einigen wenigen Tintenstrahldruckern ordnungsgemäß funktioniert. Der Ausdruck über MS-Windows erfolgt durch eine Kopie des Bildschirminhalts als Bitmap in die Zwischenablage von Windows, welches in MS- Paintbrush, MS- Powerpoint, MS- Corel Draw oder MS-Word kopiert und von dort ausgedruckt werden kann. Voraussetzung dafür ist allerdings, daß MS-Windows im Hintergrund läuft. Ist das Netzwerk zu groß, um es als Ganzes auszudrucken, kann man in NETWORK eine adäquate Vergrößerung auswählen, das Netzwerk in Teilen auf Windows übertragen und die Einzelteile z. B. unter MS- Corel Draw wieder zusammenfügen. Dieser Weg ist zwar etwas umständlich, bietet aber die Möglichkeit, aufgrund der Vorlage in Corel-Draw eine publikationsreife Zeichnung anzufertigen. Dazu müssen nur die Kanten, Kreise und Bezeichnungen des Bitmaps durch entsprechende Graphikelemente aus Corel-Draw ersetzt werden, die bzgl. Kantenstärken, Kreis- und Schriftgrößen verändert werden können.

Um einen Hinweis darauf zu geben, welche Merkmale Konflikte hervorrufen, ist eine Sektion (‚Statistics‘) enthalten, die zeigt, wie oft die Merkmale im schlechtesten Fall bei einer Auflösung des Torsos mutieren würden. Beispielsweise wird man bei RFLP- Daten stets die Schnittstelle 16517 HaeIII unter den hochmutablen Positionen finden (vgl. Forster et al. (1997)). Bei einer auffällig hohen Mutationsrate eines Merkmales bietet sich eine erneute Berechnung des Netzwerks mit einer anderen Gewichtung bzw. Entfernung des Merkmals durch Setzen des Gewichtes auf den Wert Null an.

Weiterhin wird eine Abschätzung für die Länge eines MP- Baumes geliefert. Diese ist aber nur als obere Grenze zu verstehen und keinesfalls als Wert für die exakte Länge. Der entsprechende Algorithmus berechnet zunächst innerhalb des Torsos des errechneten Netzwerkes einen kürzesten aufspannenden Baum. Danach werden Kanten, die als Endknoten Medianvektoren haben, von denen keine Endäste abzweigen, sukzessive entfernt. Zu der Länge des so entstehenden Baumes innerhalb des Torsos werden die Mutationen auf den Endästen addiert. Der so errechnete Wert dient als Schätzung der oberen Grenze für die Länge eines MP- Baumes.

Ist in dem errechneten Netzwerk kein MP- Baum enthalten, so kann der Algorithmus natürlich auch nicht die Länge eines MP- Baumes errechnen. Ebenso ist es möglich, daß der Algorithmus innerhalb des Torsos keinen kürzesten Baum der Sequenzen des Datensatzes findet, sondern nur einen kürzesten aufspannenden Baum des Torsos, in dem auch alle Medianvektoren enthalten sind, so daß die MP- Abschätzung nicht einmal für das errechnete Netzwerk optimal ist. Allgemein läßt sich aber sagen, daß die Chance für den Algorithmus, eine gute Schätzung für den kürzesten Baum innerhalb des Netzwerkes zu liefern, um so größer wird, je mehr Ausgangssequenzen des Datensatzes im Torso enthalten sind, da Kanten mit Sequenzen des Datensatzes als Endknoten bevorzugt ausgewählt werden.

 

Star Contraction (SC) Algorithmus

 

Der SC- Algorithmus kann die gleichen Datenformate verarbeiten wie der MJ- Algorithmus und ist in genau der Art und Weise implementiert, wie zuvor beschrieben.

Während des Programmlaufs werden zwei Dateien generiert. Eine Protokolldatei mit der Dateierweiterung ‘pro’, in der sämtliche Veränderungen des Datensatzes dokumentiert sind, und eine Datei mit der Erweiterung ‘sco’, die den resultierenden Datensatz in MJ- Format enthält. Die während des Programmlaufs eventuell generierten Medianvektoren werden mit ‘SCO’ zuzüglich einer Nummer benannt. Man sollte es insofern vermeiden, bereits in der Eingabedatei eine Sequenz beispielsweise mit dem Namen ‘SCO1’ zu versehen.

 

Cycle Pruning (CP) Algorithmus

 

Als Eingabeformat benötigt der CP- Algorithmus ebenso wie der Graphikalgorithmus das Ausgabeformat des MJ- Algorithmus, d. h. die Datei mit der Liste der gleichen Sequenzen, den Kanten, den Frequenzen und der Gewichtung. Nach erfolgter ‘Ausdünnung‘ des Netzwerks wird eine Ausgabedatei mit der Erweiterung ‘CYP’ generiert, die dem Ausgabeformat des MJ- Algorithmus entspricht, so daß das ‘ausgedünnte’ Netzwerk wieder mit dem Zeichenalgorithmus bearbeitet werden kann.

Fast alle Parameter sind mit Standardwerten. Einzige Ausnahme ist der Wert h . Dieser muß iterativ von einem gegebenen Startwert verringert werden, bis zum Abbruch durch den Benutzer oder bis keine weitere Ausdünnung des Netzwerks mehr möglich ist.

 

Literatur

 

Anderson S, Bankier AT, Barrell BG, de Bruijn MHL, Coulson AR, Drouin J, Eperon IC et al. (1981). Sequence and organization of the human mitochondrial genome. Nature 290: 457-465

 

Aris-Brosou S, Excoffier L (1996). The impact of population expansion and mutation rate heterogeneity on DNA sequence polymorphism. Mol Biol Evol 13: 494-504

 

Bandelt H-J, Forster P, Sykes BC, Richards MB (1995). Mitochondrial portraits of human populations using median networks. Genetics 141: 743-753

 

Bandelt H-J, Forster P, Röhl A (1999). Median joining networks for inferring intraspecific phylogenies. Mol Biol Evol 16: 37-48

 

Chen YS, Torroni A, Excoffier L, Santachiara-Benerecetti AS, Wallace DC (1995). Analysis of mtDNA variation in African populations reveals the most ancient of all human continent-specific haplogroups. Am J Hum Genet 57: 133-149

 

Excoffier L, Smouse PE, Quattro JM (1992). Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data. Genetics 131: 479-491

 

Excoffier L, Smouse PE (1994). Using allele frequencies and geographic subdivision to reconstruct gene trees within a species: Molecular variance parsimony. Genetics 136: 343-359

 

Farris JS (1970). Methods for computing Wagner trees. Syst Zool 19: 83-92

 

Fitch WM (1971). Towards defining a course of evolution: minimum change for a specific tree topology. Syst Zool 20: 406-416

 

Forster P, Harding R, Torroni A, Bandelt H-J (1996). Origin and evolution of Native American mtDNA variation: a reappraisal. Am J Hum Genet 59: 935-945

 

Forster P, Harding R, Torroni A, Bandelt H-J (1997). Reply to Bianchi and Bailliet. Am J Hum Genet 61: 245-247

 

Foulds LR, Hendy MD, Penny D (1979). A graph theoretic approach to the development of minimal phylogenetic trees. J Mol Evol 13: 127-149

 

Gusfield D (1997). Algorithms on strings, trees and sequences. Cambridge University Press, Cambridge

 

Jeffreys AJ, Wilson V, Thein SL (1985). Hypervariable ‘minisatellite’ regions in human DNA. Nature 314: 67-73

 

Jobling MA, Pandya A, Tyler-Smith C (1997). The Y chromosome in forensic analysis and paternity testing. Int J Legal Med 110: 118-124

Karp RM (1972). Reducibility among combinatorial problems. In: Complexity of computer computations (Eds. R. E. Miller und J. W. Thatcher). Plenum Press, New York: 85-103

 

Kittles RA, Perola M, Peltonen L, Bergen AW, Aragon RA, Virkkunen M, Linnoila M, Goldman D, Long JC (1998). Dual origins of Finns revealed by Y chromosome haplotype variation. Am J Hum Genet 62: 1171-1179

 

Krings M, Stone A, Schmitz RW, Krainitzki H, Stoneking M, Pääbo S (1997). Neandertal DNA sequences and the origin of modern humans. Cell 90: 19-30

 

Kruskal JB (1956). On the shortest spanning subtree of the graph and the travelling salesman problem. Proc AMS 7: 48-57

 

Lawler EL (1976). Combinatorial optimization: Networks and matroids. Holt, Rinehart and Winston, New York

 

Li WH, Graur D (1991). Fundamentals of molecular evolution. Sinauer Associates, Massachusetts

 

Meyer E, Wiegand P, Brinkmann B (1995). Phenotype differences of STRs in 7 human populations. Int J Legal Med 107: 314-322

 

Saitou N, Yamamoto F (1997). Evolution of primate AB0 blood group genes and their homologous genes. Mol Biol Evol 14: 399-411

 

Torroni A, Miller JA, Moore LG, Zamudio S, Zhuang J, Droma T, Wallace DC (1994). Mitochondrial DNA analysis in Tibet: Implications for the origin of the Tibetan population and its adaptation to high altitude. Am J Phys Anthropol 93: 189-199

 

Vigilant L, Stoneking M, Harpending H, Hawkes K, Wilson AC (1991). African populations and the evolution of human mitochondrial DNA. Science 253: 1503-1507

 

Wakeley J, Hey J (1997). Estimating ancestral population parameters. Genetics 145: 847-855

 

Watson E, Forster P, Richards M, Bandelt H-J (1997). Mitochondrial footprints of human expansions in Africa. Am J Hum Genet 61: 691-704

 

Weber JL, Wong C (1993). Mutation of human short tandem repeats. Hum Mol Genet 8: 1123-1128

 

Wills C (1995). Topiary pruning and weighting reinforce an African origin for the human mitochondrial DNA tree. Evolution 50: 977-989

 

 

Anhang B: Rolf B, Röhl A, Forster P, Brinkmann B. On the Genetic origins of the Turks

 On the Genetic Origins of the Turks

 

Burkhard Rolf1, Arne Röhl2, Peter Forster1,2, Bernd Brinkmann1

 

1Institut für Rechtsmedizin, Westfälische Wilhelms-Universität Münster, Von-Esmarch-Straße 62, D-48129 Münster;

2Mathematisches Seminar, Universität Hamburg, Bundesstraße 55, D-20146 Hamburg, Germany

 

Summary

Although language and history bear testimony to the migration of Asians to Turkey in mediaeval times, genetic evidence for this event has been elusive. Here we analyse the six Y-chromosomal short tandem repeats (STRs) DYS19, DXYS156-Y, DYS390, DYS391, DYS392, and DYS393 for 298 Turks, Germans, Japanese, and Chinese using the new median-joining (MJ) phylogenetic network method specifically designed for population DNA genealogies. The Y-chromosomal results agree with an analysis of mtDNA markers for 1210 Eurasians, indicating about 10% of east Asian genetic input into the population of Turkey.

 

Introduction

The Altaic language family extends from northeastern and central Asia to modern Turkey, forming a wedge between the Indoeuropean languages of Europe to the north and the Middle East to the south. The reason for the elongated geographic distribution to the west is known from historical sources: migrating Turkic speakers gained control of Asia minor (present Turkey), imposing their Altaic language onto the native, mainly Indoeuropean-speaking population. Given these linguistic and historical facts, it is perhaps surprising that genetic evidence for east Asian admixture into the population of Turkey is elusive. Cavalli-Sforza et al. (1994) have found that allele frequency data classify the Turks as Caucasoid, and that a more sensitive genetic tool seems to be required to detect Asian input. To this end, we have chosen Y-STR loci and mitochondrial DNA (mtDNA) variation as markers for migratory events. Both genetic systems are characterised by uniparental inheritance and a high mutation rate, affording a maximum opportunity for genetic drift after population separations. In addition, both molecules are non-recombining, permitting a phylogeographic analysis of DNA variants, which is a more delicate type of analysis than distance-based summary statistics (Richards et al. 1997).

 

Subjects and Methods

Blood samples were obtained from 80 Germans from Westphalia and surroundings (10 of these had Slavic, i.e. east German or Polish surnames: Ger3, 15, 24, 33, 48, 54, 56, 59, 71, and 75), 86 Turks from the region around Adana, 2 further Turks (Tur87 and 88) living in Germany, 72 Han Chinese from Shenyang, and 58 Japanese from the Shiga area. Primer sequences and PCR conditions for Y-chromosome typing of the STR loci DYS19, DYS390, DYS391, DYS392, and DYS393 were performed as previously described (Kayser et al. 1997). Primer sequences and amplification conditions for DXYS156 were taken from Chen et al. (1994). The DXYS156 primers amplify both a Y-specific and an X-specific fragment; however, the X- and Y-chromosome fragments are distinguishable as only allele lengths 6 to 10 are present in females and 11 to 14 in males. A single exception is a German male with 7 and 10 repeats, presumably corresponding to the X and the Y locus, respectively. A total of 226 females were typed (not shown) in the four populations to confirm the absence of allele lengths >10 on the X chromosome (cf. however Karafet et al. 1998).

Construction of the DNA genealogies of closely related sequences in a population sample usually requires phylogenetic network approaches to cope with the multitude of (almost) equally plausible trees incurred by homoplasy, whereas standard tree-building algorithms are only suitable for either less homoplasious data or interspecies data. Two different types of implemented phylogenetic network algorithms are currently available : the reduced median network algorithm (Bandelt et al. 1995) and the median-joining network algorithm (Bandelt et al. submitted). Both methods have been extensively studied on mtDNA data (Richards et al. 1996, Forster et al. 1996, Watson et al. 1997, Bandelt et al. submitted). Preliminary trials (not shown) however demonstrated that the rate of homoplasy is so high in Y-STRs that the network algorithms, applied singly on the whole data set, result in complex networks which include too many implausible trees. This observation is in agreement with others who either exclude STR data from their phylogenetic analysis (Hammer et al. 1997) or restrict the analysis to a subset of the data defined by a point mutation (Zerjal et al. 1997). To partially overcome the effect of character conflicts generated by homoplasy, we have devised a heuristic strategy for weeding out spurious homoplasious connections in a network in that we subject our data set to the following five stages of analysis. First, the pentanucleotide repeat locus DXYS156-Y was given twice the weight of other loci for both network algorithms on account of its apparent stability in major ethnic groups: over 90% of Europeans have 12 repeats at DXYS156-Y (data of Sajantila et al. 1996, Kayser et al. 1997, and this study), while about 85% of a sub-Saharan African sample (data of Sajantila et al. 1996) have 11 repeats, and 13 repeats so far have only been found in Asians (this study, excepting Ger53). Furthermore, the data of Sajantila et al. (1996) includes typing of the unique YAP insertion event, which reduces the number of possible MP trees for their African data to just four trees, all of which agree in demanding four parallel mutations at DYS19, but only one parallelism at DXYS156-Y. Second, we have extracted only those Y-types from our data set that occur in at least two individuals. The most recent mutations, which are less useful for phylogenetic reconstruction, have had little time to spread and are more likely to be represented by single individuals in a sample. In the presence of high homoplasy, there is a risk of connecting distinct clades via peripheral (recent) types when their ancestors are not present in the sample. This criterion also eliminates ancestral singletons which are rare due to age, but these are generally reconstructed as branching nodes in phylogenetic procedures. Third, we have subjected this smaller data set (195 individuals out of 298) to the reduced median network algorithm (Bandelt et al. 1995) in order to generate a solution space for which parallelisms are postulated according to character compatibility and frequency criteria to reduce the level of homoplasy. Fourth, we have subjected this output to the median-joining network algorithm with parameter e set to 0 (Bandelt et al. submitted), which successively connects subnetworks by generating only few consensus (median) sequences. Fifth, the full data set is augmented by the inferred ancestral types from the previous analysis and subjected to our combined network method (not shown). The data input was coded according to the one-step model for repeat changes, which gives a lower weight to one-repeat changes over multi-repeat changes, as inferred from pedigree studies (e.g. Brinkmann et al. submitted).

Results and Discussion

Y-chromosomal variation

For the Y-chromosomal analysis we have determined repeat numbers of the STR loci DYS19, DYS390, DYS391, DYS392, DYS393, and DXYS156-Y in Turks (n = 88), Germans (n=80), Han Chinese (n = 72), and Japanese (n = 58). These 298 individuals were found to represent 161 distinct Y-chromosomal types (or 162 if a dinucleotide deletion in Ger70 is counted). The network of the Y-types represented in at least two individuals is displayed in fig. 1. Of the 15 hypothesised ancestral nodes, eight actually exist as singletons in our sample. We have also confirmed that the full data set yields the same topology by subjecting the complete data set (including eight further southeast Europeans) to the same sort of analysis, but adding the 7 hypothesised ancestral types to the data set prior to running the two network algorithms. Only two of the 15 hypothesised nodes are superfluous in the total network (the empty nodes adjacent to Ger42 and Han40 in fig. 1). The network in the vicinity of these two nodes is therefore less robust. Although appropriate robusticity measures for trees or networks of homoplasious population data (akin to bootstrapping for species data) have not yet been developed (Bandelt et al. 1995), we can appreciate the reliability of the overall topology by the evident ethnogeographical structure. The network contains numerous reticulations, and strikingly, these nearly always involve DYS390 or DYS19, indicating that these loci mutate particularly rapidly. A count of mutation events (table 1) confirms that mutation rates between loci differ by up to an order of magnitude, suggesting that it is inadvisable to apply a uniform mutation rate to Y-STRs.

The node marked by an asterisk contains the majority state of the conservative DXYS156-Y in sub-Saharan Africans (Sajantila et al. 1996, Salem et al. 1996) and is therefore probably identical to or in the immediate vicinity of one of the major Eurasian founding Y-chromosomal types. The clusters emanating from this node are, as indicated by the colour/graphical coding, geographically specific for Caucasoids on the one hand and east Asians on the other. For example, the cluster defined by 11 repeats at DXYS156-Y is virtually absent in Europeans, but common in Asians (this study, Sajantila et al. 1996, Kayser et al. 1997, and Hammer, personal communication) and Africans (Sajantila et al. 1996, Kayser et al. 1997). Australians and highland Papuans also independently coalesce on this node (Forster et al. 1998). The vicinity of Jap11 also appears to include a major founding node for Asians, Europeans (fig. 1) as well as Papuans and Australians (Forster et al. 1998). There is, at the present genetic resolution of six STRs, a branch which does not seem to fit the general ethnic pattern. The mainly Turkish branch linked to the node marked by an asterisk contains an Asian (Han36), but no Germans. This Turkish branch may represent (a) original founders of the Turkish Y pool, which did not penetrate into or did not survive in central Europe, or (b) subsequent migration of males into our southern Turkish sample which did not affect central Europeans. With the given data, neither possibility can be dismissed. In scenario (b) two plausible historic explanations are worth consideration. The branch may conceivably be Asian, deriving from Turkic invaders. However, a historic Asian migration to Turkey may be expected to have resulted in a branch with more Asian representatives. The presence of Han36 in this branch does not necessarily lend support since Han36 is only grouped into this branch due to the fast locus DYS19 (table 1). The contact of southern Turkey to the Arabic world is also historically recorded (Rother 1971). A Near Eastern origin of the branch, whether historic or prehistoric, receives support from the fact that 50/153 Egyptians have Y types with 11 repeats at DXYS156-Y and 14 or less repeats at DYS19 (data of Salem et al. 1996).

The ethnically mixed node represented by Tur10 in fig. 1 is another apparent exception to the ethnic segregation of Y variation. However, this is simply due to insufficient resolution, and typing of the 5' variable segment of DYS389 for all the data (Rolf, unpublished) resolves the reticulation between Jap11 and Tur10 by assigning 5 repeats to Germans, and 4 repeats to all Asians descended from and including the node of Jap11. Six of the Turks (Turks 3, 10, 14, 39, 58, and 63) have this Asian motif, and one further Turk (Tur39) has the Asian-specific 13 repeats at DXYS156-Y, but 5 repeats at DYS389. At least one further "Asian" Turk is the singleton Tur64 (not included in fig. 1), which differs from Han8 by only one repeat at DYS391. The presence of the Turkish/Asian cluster descended from Jap11 in fig. 1 may again be explained by ancestral diversity, but the imperfect segregation of Turkish from Asian types may be more compatible with recent Asian immigration to Turkey. Arab immigration seems less likely in this case as 8 repeats at DXYS156-Y are not found in the Near Easterners of Salem et al. (1996). Hence, there is a minimum of 7/88 (8%) Turks who may indicate a male east Asian contribution to our southern Turkish sample.

 

mtDNA variation

Three published data sets for Turkish mtDNA exist (Richards et al. 1996, Calafell et al. 1996, Comas et al. 1996) but have not previously been applied to the question of female Turkic introgression into the population of Asia minor. The Eurasian mtDNA phylogeny has been extensively characterised in both the control region (using sequencing analysis) as well as by 14-enzyme RFLP analysis of the complete mtDNA genome, and Torroni and coworkers (Torroni et al. 1994, Torroni et al. 1996) have introduced a widely used shorthand for labeling the branches emanating from the African mtDNA founders in the Eurasian RFLP phylogeny. At the 14-enzyme RFLP resolution, the west European and east Asian branches are, unsurprisingly, almost completely disjoint phylogenetically from each other, and this is reflected in the mtDNA nomenclature (table 2): only group X (represented mainly in Amerinds) as well as groups H, U, K, T, and M* are shared between Asians and Europeans, while the other mtDNA groups are geographically specific. A proportion of 9% of the Turkish sample has east Asian mitochondrial types, and the Turkish mtDNA profile matches the Altaic-speaking Mongolians more closely than the other Asians (table 2). Due to the circumstance that up to 18 individuals (groups H, U, V, K, T, and unclassified) in the Mongolian sample are potentially indistinguishable from Caucasoid mtDNA, a figure of up to 11% of east Asian mtDNA input in the Turks may be more realistic.

 

Conclusion

This study has identified east Asian Y-chromosomal and mtDNA markers at a proportion of about 10% in the population of Turkey, confirming the conclusion of Cavalli-Sforza et al. (1994) that language replacement has occurred there with little gene replacement. Further studies using these markers can now investigate whether there are microgeographic patterns of Asian input in the Turkic-speaking populations. From a methodological point of view, it is evident that Y-STRs do allow the approximate reconstruction of DNA genealogies of a time depth of several tens of thousands of years, provided some appropriate weighting is applied. However, STR genealogies will never be as free of homoplasy as point mutation genealogies are (e.g. Hammer 1995, Jobling et al. 1997, Underhill et al. 1997) due to the high mutation rate of STRs (relative to the time depth of human evolution) coupled with locus-specific mutation rate heterogeneity, and a combination of both types of data is desirable. The advantage of the rapidly mutating STRs, apart from their forensic value, may however be the eventual quantification of absolute mutation rates in pedigree studies (Heyer et al. 1997, Kayser et al. 1997).

 

Acknowledgements

We thank M. Hammer, M. Kayser, and C. J. Kolman for information on Mongolian data, and an anonymous reviewer for valuable comments. This work was performed by A. Röhl in partial fulfilment of the requirements for the doctoral degree.

 

References

Bandelt H-J, Forster P, Sykes BC, Richards MB (1995) Mitochondrial portraits of human populations using median networks. Genetics 141:743-753

 

Bandelt H-J, Forster P, Röhl A (1999) Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol 16: 37-48

 

Calafell F, Underhill P, Tolun A, Angelicheva D, Kalaydjieva L (1996) From Asia to Europe: mitochondrial DNA sequence variability in Bulgarians and Turks. Ann Hum Genet 60:35-40

 

Cavalli-Sforza LL, Menozzi P, Piazza A (1994) The History and Geography of Human Genes. Princeton University Press, New Jersey, p 245

 

Chen H, Lowther W, Avramopoulos D, Antonarakis SE (1994) Homologous loci DXYS156X and DXYS156Y contain a polymorphic pentanucleotide repeat (TAAAA)n and map to human X and Y chromosomes. Human Mutation 4:208-211

 

Comas D, Calafell F, Mateu E, Pérez-Lezaun A, Bertranpetit J (1996) Geographic variation in human mitochondrial control region sequence: the population history of Turkey and its relationship to the European populations. Mol Biol Evol 13:1067-1077

 

Forster P, Harding R, Torroni A, Bandelt H-J (1996) Origin and evolution of Native American mtDNA variation: a reappraisal. Am J Hum Genet 59:935-945

 

Forster P, Kayser M, Meyer E, Roewer L, Pfeiffer H, Benkmann H, Brinkmann B (1998) Phylogenetic resolution of complex mutational features at Y-STR DYS390 in aboriginal Australians and Papuans. Mol Biol Evol (in press)

 

Hammer MF (1995) A recent common ancestry for human Y chromosomes. Nature 378:376-378

 

Hammer MF, Spurdle AB, Karafet T, Bonner MR, Wood ET, Novelletto A, Malaspina P et al (1997) The geographic distribution of human Y chromosome variation. Genetics 145:787-805

 

Heyer E, Puymirat J, Dieltjes P, Bakker E, de Knijff P (1997) Estimating Y chromosome specific microsatellite mutation frequencies using deep rooting pedigrees. Hum Mol Genet 6:799-803

 

Horai S et al (1996) mtDNA polymorphism in east Asian populations, with special reference to the peopling of Japan. Am J Hum Genet 59:579-590

 

Jobling MA, Pandya A, Tyler-Smith C (1997) The Y chromosome in forensic analysis and paternity testing. Int J Legal Med 110:118-124

 

Karafet T, de Knijff P, Wood E, Ragland J, Clark A, Hammer MF (1998) Different patterns of variation at the X- and Y-linked microsatellite loci DXYS156X and DXYS156Y in human populations. Hum Biol (in press)

 

Kayser M, Caglià A, Corach D, Fretwell N, Gehrig C, Graziosi G, Heidorn F et al (1997) Evaluation of Y-chromosomal STRs: a multicenter study. Int J Legal Med 110:125-133 and 141-149

 

de Knijff P, Kayser M, Caglià A, Corach D, Fretwell N, Gehrig C, Graziosi G et al (1997) Chromosome Y microsatellites: population genetic and evolutionary aspects. Int J Legal Med 110:134-140

 

Kolman CJ, Sambuughin N, Bermingham E (1996) Mitochondrial DNA analysis of Mongolian populations and implications for the origin of New World founders. Genetics 142:1321-1334

 

Richards, MB, Macaulay V, Sykes BC, Pettitt P, Forster P, Bandelt HJ (1997) Reply to Cavalli-Sforza and Minch. Am J Hum Genet 61:251-253

 

Richards MB, Côrte-Real H, Forster P, Macaulay V, Wilkinson-Herbots H, Demaine A, Papiha S et al (1996) Paleolithic and neolithic lineages in the European mitochondrial gene pool. Am J Hum Genet 59:185-203

 

Rother L (1971) Die Städte der Çukurova: Adana-Mersin-Tarsus. Geographisches Institut der Universität Tübingen, Tübingen

 

Ruhlen M (1987) A Guide to the World's Languages. Stanford University Press, Stanford, California

 

Sajantila A, Salem A-H, Savolainen P, Bauer K, Gierig C, Pääbo S (1996) Paternal and maternal DNA lineages reveal a bottleneck in the founding of the Finnish population. Proc Natl Acad Sci USA 93:12035-12039

 

Salem A-H, Badr FM, Gaballah MF, Pääbo S (1996) The genetics of traditional living: Y-chromosomal and mitochondrial lineages in the Sinai peninsula. Am J Hum Genet 59:741-743

 

Torroni A, Miller JA, Moore LG, Zamudio S, Zhuang J, Droma T, Wallace DC (1994) Mitochondrial DNA analysis in Tibet: implications for the origin of the Tibetan population and its adaption to high altitude. Am J Phys Anthropol 93:189-199

 

Torroni A, Huoponen K, Francalacci P, Petrozzi M, Morelli L, Scozzari R, Obinu D et al (1996) Classification of European mtDNAs from an analysis of three European populations. Genetics 144:1835-1850

 

Underhill PA, Jin L, Lin AA, Quasim Mehdi S, Jenkins T, Vollrath D, Davis RW et al (1997) Detection of numerous Y chromosome biallelic polymorphisms by denaturing high-performance liquid chromatography. Genome Res 7:996-1005

 

Watson E, Forster P, Richards MB, Bandelt H-J (1997) Mitochondrial footprints of human expansions in Africa. Am J Hum Genet 61:691-704

 

Zerjal T, Dashnyam B, Pandya A, Kayser M, Roewer L, Santos FR, Schiefenhövel W et al (1997) Genetic relationships of Asians and northern Europeans, revealed by Y-chromosomal DNA analysis. Am J Hum Genet 60:1174-1183

 

Table 1

Estimated number of mutations at Y STR loci

 

minimum number of mutations

maximum number of mutations

DYS19

15

28

DXYS156-Y

3

(3)

DYS390

13

28

DYS391

7

12

DYS392

7

15

DYS393

7

10

 

Note: The estimated number of mutations for each locus was determined from the phylogenetic network in fig. 1 (58 Y-types). Minimum and maximum mutations were counted for a given locus by resolving reticulations respectively in favour of or against a given locus. The maximum number of mutations for DXYS156-Y is not comparable since this locus is weighted for the phylogeny (see Subjects and Methods).

 

Table 2

mtDNA groups in Europeans, Turks, and Asians

 

 

Asian mtDNA groups

European mtDNA groups

African

unclass.

 

n

A

B

C

D

E

F

M*

H,U,V

I,W,X

J

K

T

L1,L2

 

W. Europe

891

0

0

0

0

0

0

1

596

44

105

61

75

7

2

Turkey

96

1

0

3

4

0

0

1

50

10

15

3

8

1

0

Mongolia

103

4

10

15

25

2

9

20

9

0

0

1

1

0

7

Tibet

54

6

3

2

9

4

9

21

0

0

0

0

0

0

0

Taiwan

66

5

13

2

7

2

15

16

1

0

0

0

0

0

5

 

Notes: European data from Richards et al. (1996) and Torroni et al. (1996); Turkish data from Richards et al. (1996), Calafell et al. (1996) and Comas et al. (1996); Mongolian data from Kolman et al. (1996); Tibetan data from Torroni et al. (1994); Taiwanese (Han) data from Horai et al. (1996). mtDNA group nomenclature according to Torroni et al. (1996) and Bandelt et al. (submitted). mtDNA lineages M*, C, D, E together constitute the supergroup M.

 

 

 

Figure 1 Phylogenetic network of Eurasian Y-STR types. Only types that occur at least twice in the data set (195 individuals out of the total data set of 298) are displayed. Types shared between populations are labelled by only one individual. Circle area is proportional to the number of individuals, and ethnic origins are depicted with colours as indicated. The node marked by an asterisk represents corresponds to the following sequence (amplicon length indicated in parentheses): DYS390: 24 repeats (215bp), DYS19: 15 repeats (194bp), DYS392: 11 repeats (248bp), DYS393: 13 repeats (124bp), DYS391: 10 repeats (279bp), DXYS156-Y: 11 repeats (160bp). Arrows point to the direction of a repeat gain.

 

 

Zusammenfassung

 

Bei der phylogenetischen Analyse nicht-rekombinanter intraspezifischer DNS-Daten sind die (minimalen) Steinerbäume im Sequenzraum (bzgl. der Hammingdistanz) zu bestimmen. Für drei Sequenzen läßt sich das Steinerproblem trivial mittels der quasimedianen Operation lösen, die einen Medianvektor gemäß Majoritäts- und Prioritätsregel als Steinerpunkt bestimmt. Durch Iterieren dieser Operation erhält man eine Sequenzmenge, die, als (gewichteter) Graph ("quasimedianes Netzwerk") organisiert, einen universellen Lösungsraum für das Steinerproblem darstellt. Da für typische mitochondriale DNS-Datensätze angesichts zahlreicher Rück- und Parallelmutationen das volle quasimediane Netzwerk jedoch zu komplex ist, stellt sich das Problem, Heuristiken zu entwickeln, die nur lokal einfach zu bestimmende Teile des quasimedianen Netzwerks zu einem weniger komplexen Netzwerk zusammenfügen und dieses schließlich zu einem fast baumhaften Netzwerk ausdünnen, das aus mutmaßlichen Steinerlösungen besteht, die im Hinblick auf zusätzliche Kriterien als vielversprechende Kandidaten für die gesuchte wahre Phylogenie angesehen werden. Insbesondere werden die beiden folgenden Heuristiken entwickelt:

  • Star Contraction zerlegt den Datensatz in "sternförmige" Strukturen.
  • Cycle Pruning entfernt sukzessive Kanten aus einem Teilnetzwerk des quasimedianen Netzwerks unter Berücksichtigung der Häufigkeit gesammelter Sequenzen im gegebenen Datensatz.

Diese beiden Algorithmen werden in Verbindung mit zwei bereits publizierten Netzwerkverfahren in einer Simulationsstudie getestet.

Im Anhang A wird das Programmpaket NETWORK 2.0 beschrieben, und im Anhang B wird eine Fallstudie zur Anwendung von Netzwerkverfahren auf Y-chromosomale STRs (short tandem repeats) wiedergegeben.

 

 

* Lebenslauf (Curriculum vitae)

 

Allgemeines

Name

Arne Röhl

Geburtsdatum

08.04.1969

Familienstand

Ledig

Eltern

Dr. jur. Ernst Adolf Röhl

Eveline Röhl

Ausbildung

Grundschule

1975-1979 in Prisdorf bei Pinneberg

Gymnasium

1979-1988 an der Johannes Brahms Schule in Pinneberg

Studium der Mathematik

1989-1996 am Mathematischen Seminar der Universität Hamburg

Promotion

1996-1999 bei Prof. Dr. H. J. Bandelt, Universität Hamburg

Beruf

1988-1989

Wehrdienst in Wentorf bei Hamburg

Februar 1997 – Juni 1997

 

Wissenschaftlicher Mitarbeiter des Mathematischen Seminars der Universität Hamburg

Juli 1997 – September 1998

Mathematischer Mitarbeiter bei der KRAVAG-Versicherung, Hamburg

Seit November 1998

Consultant bei der Unternehmensberatung PricewaterhouseCoopers, Geschäftsstelle Hamburg

 

  

 

 

 

Erläuterungen:

 

Ein Netzwerk ist definiert als die Vereinigung endlich vieler ungerichteter Bäume der gleichen Knoten und kann dadurch im Gegensatz zu einem Baum Kreise enthalten.

Zurück

 

Zwei binäre Merkmale mit den Zuständen 0 und 1 sind genau dann inkompatibel, wenn es vier Sequenzen gibt, die die Kombinationen 00, 01, 10 und 11 bzgl. der beiden Merkmale enthalten, d. h. daß es nicht möglich ist, die Merkmale in einem Baum zu repräsentieren in dem jeweils genau eine Mutation vorkommt.

Zurück

 

Ein Kreis liegt vor, wenn es einen Kantenzug (e1,...,en) paarweise verschiedener Kanten ei und Punkte v0,...,vn mit ei=vi-1vi, i=1,...,n, innerhalb des Graphen gibt, für die gilt v0=vn.

Zurück

 

Ein dichotomer Baum zeichnet sich dadurch aus, daß jeder innere Knoten genau mit drei Kanten inzidiert, bzw. sich der beliebig gewurzelte Baum an jeder Gabelung in zwei Äste aufspaltet.

Zurück

 

Eine Kante zwischen den Sequenzen x und y ist genau dann ‘e -zulässig’, wenn es keinen aus Sequenzen des Datensatzes bestehenden Weg x = u0, u1,..., uk = y zur Verbindung von x und y gibt, so daß gilt: d(ui ,ui+1 ) < d(x , y) - e für alle i=0,...,k-1.

Zurück

 

Der Torso eines Netzwerkes besteht aus allen Verbindungen, die in mindestens einem Kreis enthalten sind, sowie allen Verbindungen zwischen Kreisen.

Zurück