Informație

8.7: Anexă - Algoritmul de tăiere al lui Felsenstein - Biologie


Algoritmul de tăiere al lui Felsenstein (1973) este un exemplu de programare dinamică, un tip de algoritm care are multe aplicații în biologia comparativă. În programarea dinamică, descompunem o problemă complexă într-o serie de pași mai simpli care au o structură imbricată. Acest lucru ne permite să refolosim calculele într-un mod eficient și accelerăm timpul necesar efectuării calculelor.

Cel mai bun mod de a ilustra algoritmul lui Felsenstein este printr-un exemplu, care este prezentat în panourile de mai jos. Încercăm să calculăm probabilitatea unui caracter cu trei stări pe un arbore filogenetic care include șase specii.

Figura 8.4A. Fiecare vârf și nod intern din arbore are trei casete, care vor conține probabilitățile pentru cele trei stări de caracter în acel moment al arborelui. Prima casetă reprezintă o stare 0, a doua stare 1 și a treia stare 2. Imaginea autorului, poate fi reutilizată sub o licență CC-BY-4.0.

  1. Primul pas în algoritm este să completați probabilitățile pentru sfaturi. În acest caz, cunoaștem stările de la vârfurile copacului. Matematic, afirmăm că cunoaștem cu precizie stările de caracter la vârfuri; probabilitatea ca acea specie să aibă starea pe care o observăm este 1, iar toate celelalte stări au probabilitatea zero:
  1. Apoi, identificăm un nod în care toți descendenții săi imediați sunt sfaturi. Va exista întotdeauna cel puțin un astfel de nod; de multe ori, vor fi mai multe, caz în care vom alege în mod arbitrar una. Pentru acest exemplu, vom alege nodul care este cel mai recent strămoș comun al speciilor A și B, etichetat ca nod 1 în Figura 8.2B.
  2. Apoi folosim ecuația 7.6 pentru a calcula probabilitatea condițională pentru fiecare stare de caracter pentru subarborele care include nodul pe care l-am ales în pasul 2 și descendenții săi de vârf. Pentru fiecare stare de caracter, probabilitatea condiționată este probabilitatea, date fiind datele și modelul, de a obține stările de caracter tip dacă începeți cu acea stare de caracter la rădăcină. Cu alte cuvinte, ținem evidența probabilității pentru părțile descendente ale copacului, inclusiv datele noastre, dacă nodul pe care îl luăm în considerare avea fiecare dintre stările de caracter posibile. Acest calcul este:

$$ L_P (i) = ( sum limits_ {x in k} Pr (x | i, t_L) L_L (x)) cdot ( sum limits_ {x in k} Pr (x | i, t_R) L_R (x)) label {8.7} $$

Unde eu și X sunt ambii indici pentru k stări de caracter, cu sume luate în toate stările posibile la vârfurile ramurii (X) și termenii calculați pentru fiecare stare posibilă la nod (eu). Cele două piese ale ecuației sunt descendentul stâng și dreapta al nodului de interes. Ramurile pot fi atribuite în mod arbitrar la stânga sau la dreapta fără a afecta rezultatul final, iar abordarea funcționează și pentru politomii (dar ecuația este ușor diferită). Mai mult, fiecare dintre aceste două piese în sine are două părți: probabilitatea de a începe și de a se termina cu fiecare stare de-a lungul celor două ramuri fiind luate în considerare și probabilitățile condiționale actuale care intră în ecuație la vârfurile subarborelui (LL(X) și LR(X)). Lungimile ramurilor sunt indicate ca tL și tR pentru stânga și respectiv pentru dreapta.

Ne putem gândi la probabilitatea de „curgere” pe ramurile arborelui, iar probabilitățile condiționate pentru ramurile stânga și dreapta se combină prin multiplicare la fiecare nod, generând probabilitatea condițională pentru nodul părinte pentru fiecare stare de caracter (LP(eu)).

Luați în considerare subarborele care duce la speciile A și B în exemplul dat. Cele două stări de caracter tip sunt 0 (pentru speciile A) și 1 (pentru speciile B). Putem calcula probabilitatea condițională pentru starea de caracter 0 la nodul 1 ca:

$$ L_P (0) = left ( sum limits_ {x in k} Pr (x | 0, t_L = 1.0) L_L (x) right) cdot left ( sum limits_ {x in k} Pr (x | 0, t_R = 1.0) L_R (x) right) label {8.8} $$

Apoi, putem calcula termenii de probabilitate din matricea de probabilitate P. În acest caz tL = tR = 1.0, deci atât pentru ramura stângă, cât și pentru dreapta:

$$ mathbf {Q} t = begin {bmatrix} -2 & 1 & 1 1 & -2 & 1 1 & 1 & -2 end {bmatrix} cdot 1.0 = begin { bmatrix} -2 & 1 & 1 1 & -2 & 1 1 & 1 & -2 end {bmatrix} label {8.9} $$

Astfel încât:

$$ mathbf {P} = e ^ {Qt} = begin {bmatrix} 0.37 & 0.32 & 0.32 0.32 & 0.37 & 0.32 0.32 & 0.32 & 0.37 end {bmatrix} label {8.10} $$

Acum observați că, din moment ce starea caracterului stâng este cunoscută ca fiind zero, LL(0) = 1 și LL(1)=LL(2) = 0. În mod similar, starea potrivită este una, deci LR(1) = 1 și LR(0)=LR(2)=0.

Acum putem completa cele două părți ale ecuației 8.2:

$$ sum limits_ {x in k} Pr (x | 0, t_L = 1.0) L_L (x) = 0.37 cdot 1 + 0.32 cdot 0 + 0.32 cdot 0 = 0.37 label {8.11} $$

și:

$$ sum limits_ {x in k} Pr (x | 0, t_R = 1.0) L_R (x) = 0.37 cdot 0 + 0.32 cdot 1 + 0.32 cdot 0 = 0.32 label {8.11B} $ $

Asa de:

LP(0)=0.37 ⋅ 0.32 = 0.12.

Aceasta înseamnă că, sub model, dacă starea de la nodul 1 ar fi 0, am avea o probabilitate de 0,12 pentru această mică secțiune a arborelui. Putem folosi o abordare similară pentru a constata că:

(ec. 8.13)

LP(1)=0.32 ⋅ 0.37 = 0.12.

LP(2)=0.32 ⋅ 0.32 = 0.10.

Acum avem probabilitatea pentru toate cele trei stări ancestrale posibile. Aceste numere pot fi introduse în casetele corespunzătoare:

  1. Apoi repetăm ​​calculul de mai sus pentru fiecare nod din copac. Pentru nodurile 3-5, nu toate LL(X) și LR(X) termenii sunt zero; valorile lor pot fi citite din casetele de pe copac. Rezultatul tuturor acestor calcule:
  1. Acum putem calcula probabilitatea pe tot arborele folosind probabilitățile condiționate pentru cele trei stări de la rădăcina arborelui.

$$ L = sum limits_ {x in k} pi_x L_ {root} (x) label {8.14} $$

Unde πX este probabilitatea anterioară a stării de caracter la rădăcina arborelui. Pentru acest exemplu, vom lua aceste probabilități anterioare ca fiind uniforme, egale pentru fiecare stat (πX = 1/k = 1/3). Probabilitatea pentru exemplul nostru, atunci, este:

[L = 1/3 ⋅ 0.00150 + 1/3 ⋅ 0.00151 + 1/3 ⋅ 0.00150 = 0.00150 etichetă {8.15} ]

Rețineți că, dacă încercați acest exemplu într-un alt pachet software, cum ar fi GEIGER sau PAUP *, software-ul va calcula o probabilitate ln de -6,5, care este exact jurnalul natural al valorii calculate aici.


O problemă de reconstrucție pentru o clasă de rețele filogenetice cu transferuri laterale de gene

Transferurile genetice laterale sau orizontale sunt un tip de evenimente evolutive asimetrice în care materialul genetic este transferat de la o specie la alta. În această lucrare luăm în considerare Rețele LGT, un model general de rețele filogenetice cu transferuri de gene laterale care constă, aproximativ, dintr-un principal copac înrădăcinat cu frunzele etichetate pe un set de taxoni și un set de extra secundar arcuri între noduri din acest copac reprezentând transferuri de gene laterale. O rețea LGT dă naștere într-un mod natural unui subarborele filogenetic principal și un set de subarburi filogenetici secundari, care, aproximativ, reprezintă, respectiv, linia principală de evoluție a majorității genelor și liniile secundare de evoluție prin transferuri laterale de gene.

Rezultate

Introducem un set de condiții simple pe o rețea LGT care garantează că subarborii filogenetici principali și secundari sunt diferiți în perechi și că acești subarburi determină, până la izomorfism, rețeaua LGT. Apoi dăm un algoritm care, având în vedere un set de copaci filogenetici diferiți în perechi (T_0, T_1, ldots, T_k ) pe același set de taxoni, produce, atunci când există, rețeaua LGT care îndeplinește aceste condiții și astfel încât arborele său filogenetic principal este (T_0 ), iar arborii filogenetici secundari sunt (T_1, ldots, T_k ).


Abstract

Modelele de coevolutie moleculara pot dezvălui constrângeri structurale și funcționale în interiorul sau printre moleculele organice. Aceste tipare sunt mai bine înțelese atunci când se ia în considerare procesul evolutiv de bază, care ne permite să separăm semnalul evoluției dependente a siturilor (coevolutie) de efectele strămoșilor comune ale genelor. În schimb, ignorarea evoluției dependente a siturilor atunci când se studiază istoria genelor are un impact negativ asupra acurateței arborilor filogenetici deduși. Deși coevoluția moleculară și istoria filogenetică sunt interdependente, analizele celor două procese sunt efectuate separat, o alegere dictată de comoditatea de calcul, dar în detrimentul acurateței. Prezentăm o metodă bayesiană și un software asociat pentru a deduce câte și ce site-uri ale unei aliniere evoluează în funcție de un proces evolutiv independent sau dependent de perechi și pentru a estima simultan relațiile filogenetice între secvențe. Ne validăm metoda pe seturi de date sintetice și ne provocăm predicțiile coevoluției pe molecula de ARNr 16S comparându-le cu structura sa moleculară cunoscută. În cele din urmă, evaluăm acuratețea arborilor filogenetici dedusă sub presupunerea independenței între site-uri folosind seturi de date sintetice, molecula de ARNr 16S și 10 aliniamente suplimentare ale genelor eucariote care codifică proteinele. Rezultatele noastre demonstrează că deducerea arborilor filogenetici, în timp ce se ține cont de evoluția dependentă a sitului, are un impact semnificativ asupra estimărilor filogeniei și asupra procesului evolutiv.

Coevoluția moleculară este procesul evolutiv prin care interacțiunile dintre siturile îndepărtate ale unei molecule sau locurile diferitelor molecule sunt menținute astfel încât să păstreze proprietăți funcționale sau structurale avantajoase. De exemplu, fragmentele coevolutive din cadrul secvențelor de proteine ​​sunt implicate în constrângeri de pliere și informative privind intermediarii de pliere, asamblarea peptidelor sau mutații cheie cu roluri cunoscute în bolile genetice (1, 2). Disponibilitatea în continuă creștere a secvențelor moleculare (nucleotide și aminoacizi) ne oferă o cantitate fără precedent de date care dețin un potențial puternic de a dezvălui gene și regiuni genetice care evoluează sub un proces constrâns (3, 4).

Există mai multe metode pentru a deduce coevoluția doar din datele secvenței (bazate pe modele de potrivire între site-uri, așa cum sunt revizuite în referințele 5 și 6). Cu toate acestea, aceste metode nu exploatează o componentă cheie în modelarea proceselor evolutive subiacente: arborele filogenetic care descrie relațiile dintre secvențele moleculare. Incorporarea semnalului filogenetic în analiza coevoluției este crucială, deoarece ne permite să facem distincția între modele cu adevărat coevolutive și modele similare induse de istoria comună a secvențelor (5, 7). În acest scop, au fost dezvoltate mai multe metode pentru a deduce coevoluția în timp ce contabilizează relațiile filogenetice (7), dar doar câteva dintre acestea modelează în mod explicit procesul de coevolutie de-a lungul unui anumit arbore filogenetic (8 ⇓ ⇓ –11).

Toate metodele conștiente de filogenie pentru a detecta coevolutia se bazează pe presupunerea că relațiile filogenetice dintre secvențe sunt cunoscute și pot fi tratate ca „date observate”. De obicei, arborii filogenetici sunt înșiși deduși din datele moleculare (12), dar inferența lor se bazează pe o presupunere fundamentală că fiecare sit evoluează independent de toate celelalte (13). Această presupunere, care este în mod evident încălcată în prezența coevoluției, are beneficii în ceea ce privește tractabilitatea calculațională, deoarece probabilitatea unei alinieri date unui copac filogenetic este produsul probabilității individuale a fiecărui sit. S-a demonstrat că această simplificare a mecanismului evolutiv în prezența unor site-uri nedependente scade precizia arborilor filogenetici deduși (14, 15). Cu toate acestea, seturile de date cu constrângeri funcționale sau structurale puternice sunt adesea analizate în cadrul filogenetic care își asumă independența între site-uri. De exemplu, subunitatea ribozomală mică (16S) este frecvent utilizată pentru a estima cele mai vechi relații evolutive dintre descendențele majore ale arborelui vieții (16, 17), neglijând numeroasele sale constrângeri structurale și dovezile coevoluției (18).

Prezența modelelor coevolutive în numeroase secvențe de nucleotide și aminoacizi se extinde cu mult dincolo de gena 16S și este susținută de un număr mare de dovezi (9, 19). Ignorarea interdependențelor dintre istoria filogenetică a secvențelor și procesele constrânse care guvernează evoluția nucleotidelor sau aminoacizilor poate împiedica grav capacitatea noastră de a deduce arborii filogenetici corecți (14, 15) și de a detecta cu precizie coevolutia (5, 7). Cu toate acestea, inferența acestor procese interdependente este încă efectuată separat pentru comoditate matematică, în detrimentul asumării independenței între site-uri care a fost inițial descrisă ca „nu neapărat biologic valabilă” de Felsenstein în 1983, în zorii filogeneticii moleculare bazate pe probabilitate ( 13). Pentru a aborda această problemă, vă prezentăm un cadru bayesian pentru a analiza o aliniere a nucleotidelor și a estima în comun (eu) numărul de perechi de site-uri care evoluează împreună și poziția lor în secvență, diferențându-se astfel de un model de bază de evoluție independentă, (ii) parametrii modelelor independente și dependente de substituție și (iii) arborele filogenetic care descrie relațiile dintre secvențe. Metoda noastră se numește CoevRJ, iar software-ul care o implementează este disponibil la https://bitbucket.org/XavMeyer/coevrj.

Am evaluat performanța CoevRJ în reconstrucția copacilor filogenetici și deducerea coevoluției pe baza unei game extinse de seturi de date simulate, o aliniere a ARNr-ului 16S foarte coevolut și a 10 seturi de date empirice eucariote ale genelor care codifică proteinele. Arătăm că CoevRJ oferă o identificare exactă a site-urilor care evoluează, astfel cum sunt validate de date empirice și simulate. Evaluăm efectele coevoluției asupra estimărilor filogenetice prin compararea rezultatelor noastre cu cele obținute în ipoteza unei evoluții independente și demonstrăm importanța contabilității dependenței între site-uri atunci când se deduce arborii filogenetici pe seturile de date supuse coevoluției.


Metode bazate pe distanță

După citirea în aliniere putem construi un prim arbore cu metode bazate pe distanță. Funcția dist.dna din pachetul maimuță calculează distanțele pentru multe modele de substituție ADN. Pentru a utiliza funcția dist.dna trebuie să transformăm datele în clasa DNAbin. Pentru aminoacizi funcția dist.ml oferă modele comune de substituție (de exemplu „WAG”, „JTT”, „LG”, „Dayhoff”, „cpREV”, „mtmam”, „mtArt”, „MtZoa” sau „mtREV24” ).

După construirea unei matrice de distanță, reconstituim un copac înrădăcinat cu UPGMA și, alternativ, un copac neînrădăcinat folosind Neighbor Joining (Saitou și Nei 1987) (Studier și Keppler 1988). Mai multe metode de distanță precum fastme sunt disponibile în maimuţă pachet.

Putem parcela arborii treeUPGMA și treeNJ cu comenzile:

Metodele bazate pe distanță sunt foarte rapide și vom folosi arborele UPGMA și NJ ca arbori de pornire pentru analize de maximă parcimonie și maximă probabilitate.


Concluzie

Secvențierea unui număr mare de genomi întregi a făcut posibilă studierea modelelor de câștig și pierdere a genelor pe o scară enormă. Chiar dacă metodele de deducere a câștigului și pierderii au existat de aproape 30 de ani, numai odată cu analiza genomului întreg efectul unor mici prejudecăți a devenit clar. Efectul net al prejudecății în metodele de reconciliere a copacilor demonstrat aici este că numărul pierderilor genetice deduse nu ar trebui luat la valoarea nominală, iar numărul de duplicări deduse ar trebui analizat cu grijă.

Cum am putea depăși părtinirea metodelor actuale de reconciliere? Algoritmii care oferă o estimare a suportului statistic pentru fiecare câștig sau pierdere dedusă (de exemplu, [12, 30]) sau care iau în considerare lungimea ramurilor de copaci ale speciilor pot oferi ambele îmbunătățiri metodelor actuale. Această din urmă posibilitate oferă o modalitate de a îmbunătăți inferențele, deoarece ramurile copacilor de specii scurte sunt cele care sunt cel mai probabil să conducă la arborii genetici deduși în mod greșit, așa cum se vede în cazul sortării incomplete a descendenței. De asemenea, este important să rețineți că majoritatea prejudecăților descrise aici apar atunci când există un număr egal de gene între taxoni - dacă există un număr inegal de gene, atunci duplicările și pierderile pot fi deduse din informațiile de prezență / absență. În acest caz, problema se reduce pur și simplu la una referitoare la evoluția numărului de copii, care este exact abordarea metodei de probabilitate menționată mai sus. Dar această metodă nu este lipsită de propriile părtiniri [3].

În cele din urmă, părtinirile în metodele de reconciliere a arborilor pun la îndoială lucrările anterioare în evoluția genomului vertebratelor. Cu toate acestea, rezultatele prezentate aici nu pot respinge rezultatele anterioare, deoarece un exces de duplicate antice și pierderile recente de gene au fost deduse chiar și atunci când se utilizează limite de bootstrap de 100%. Descoperirea adevăratului model al evoluției genomului vertebratelor va necesita rezultate de simulare, arbori genetici mai buni, mai multe date sau o combinație a tuturor celor trei.


Concluzie

Vă prezentăm MowgliNNI metoda euristică bazată pe rearanjări NNI ale părților incerte ale arborilor genetici pentru a rezolva o problemă de optimizare a parsimoniilor pentru reconcilieri care contează duplicări (), pierderi () și transferuri (). Arătăm dovezi experimentale că reconcilierile calculate sub criteriul parsimonii pot corecta în mod eficient părți eronate ale arborilor genetici deduse din datele secvenței. Pe date simulate, MowgliNNI propune adesea o nouă topologie a arborelui genetic care este mai aproape de cea corectă și care duce, de asemenea, la mai bune și la predicții. Mai mult decât atât, numărul evenimentelor și numărul celor mai parțiale reconcilieri prezise de MowgliNNI sunt semnificativ mai mici decât cele obținute fără a pune sub semnul întrebării topologia arborelui genetic. Acest lucru este confirmat pe date reale. Un punct critic pentru metodele de parsimonie este alegerea costurilor respective pentru evenimentele evolutive considerate. Arătăm aici că MowgliNNIPerformanța este ușor modificată la schimbarea costurilor acordate evenimentelor individuale (, și), adică metoda este robustă pentru a specifica greșit costurile.


Estimarea suportului

Dat T X și T y într-o clasă de echivalență, un join T X y pe perechea comandată (T X,T y) este un FST dacă este susținut de cel puțin o fracțiune f dintre arborii de intrare (adică, dacă este cel puțin o fracțiune f dintre arborii de intrare au T X y ca subarbore). Rețineți că orice astfel de T X y este susținut doar de acei copaci care îi susțin pe amândoi T X și T y de asemenea. Motivat de aceasta, pentru fiecare FST T X menținem un lista de asistență, notat cu T X.supList, care conține toți copacii din colecția de intrări care sunt suportați T X. Pentru a estima dacă aderarea la (T X,T y) rezultă un FST, aplicăm teorema 1 numai pe copacii din T X.supList∩T y.supList. Stocăm lista de asistență ca o bitmap [36] pentru utilizarea eficientă a memoriei și calculul rapid al intersecției listelor de asistență.


8.7: Anexă - Algoritmul de tăiere al lui Felsenstein - Biologie

Partea I: ANALIZA SECVENȚEI.

1. Alinierea globală a secvențelor în perechi.

1.1 Alinierea și evoluția.

1.3 O schemă de notare pentru model.

1.4 Găsirea alinierilor cu cel mai mare scor cu programarea dinamică.

1.4.1 Determinați Hi, j.

1.4.3 Găsirea alinierilor care dau cel mai mare scor.

1.6 Lacune de scor: penalizări de lacune.

1.7 Programare dinamică pentru pedeapsa generală a decalajului.

1.8 Programare dinamică pentru penalizarea decalajului afin.

1.9 Scorul de aliniere și distanța de secvență.

2 Aliniere locală pereche și căutare în baza de date.

2.1 Operațiunea de bază: compararea a două secvențe.

2.3.2 Găsirea celor mai bune aliniamente locale.

2.3.4 Matrici de scor și penalizări gap.

2.4.2 Preprocesați interogarea: faceți lista de cuvinte.

2.4.3 Scanarea secvențelor bazei de date.

3. Analiza statistică.

3.1 Testarea ipotezei pentru omologia secvenței.

3.1.1 Generarea aleatorie de secvențe.

3.1.2 Utilizarea Z valori pentru estimarea semnificației statistice.


Metode

Filogenii ontogenetice

Modelăm mitocondriile din țesuturile prelevate de la unul sau mai mulți indivizi înrudiți ca un grup de populații legate de o filogenie ontogenetică. De-a lungul fiecărei ramuri a filogeniei ontogenetice, frecvențele heteroplasmei din interiorul unor țesuturi ancestrale se schimbă datorită acțiunii derivării genetice și a mutației. Presupunem că este dată forma filogeniei ontogenetice.

Modelul nostru de filogenie ontogenetică diferă în câteva moduri importante de cadrul tipic de probabilitate a populației-filogenetic. În modelul genetic tipic al populației, fiecare ramură este considerată a fi o perioadă independentă a istoriei evolutive și, prin urmare, se află sub controlul unui parametru independent. Spre deosebire de aceasta, permitem unui singur parametru să determine dinamica pe mai multe părți ale filogeniei, deoarece un singur proces de dezvoltare poate acționa la mai mulți indivizi înrudiți, iar acest proces de dezvoltare poate fi presupus că acționează în mod similar la fiecare individ. Mai mult, deși se presupune de obicei că fiecare locus a fost transmis printr-o singură filogenie și, prin urmare, a fost supus acelorași forțe genetice ale populației, permitem efectelor derivei genetice și ale mutației să depindă de vârsta indivizilor eșantionați. În special, pentru anumite procese ontogenetice, modelăm rata acumulării derivei genetice și a mutației odată cu vârsta. Acest lucru este motivat de observațiile anterioare că variantele heteroplasmatice se separă și se acumulează în timp în țesuturile somatice (Sondheimer și colab. 2011 Rebolledo-Jaramillo și colab. 2014 Li și colab. 2015) și în cadrul liniei germinale (Rebolledo-Jaramillo și colab. 2014 Li și colab. 2016 Wachsmuth și colab. 2016). În cele din urmă, în modelul tipic de populație-filogenetic, fiecare ramură a filogeniei reprezintă o singură perioadă din istoria evoluției și este modelată de un singur parametru. Deoarece multiple procese ontogenetice de interes pot apărea de-a lungul unei singure ramuri a filogeniei ontogenetice, permitem ramificărilor filogeniei ontogenetice să fie împărțite în multiple procese ontogenetice distincte, controlate de parametri independenți. Figura 1 demonstrează aceste caracteristici cu o filogenie ontogenetică care reprezintă relațiile dintre două țesuturi prelevate atât la mamă, cât și la descendenții ei.

Filogenia ontogenetică pentru țesuturile prelevate în duo-mamă-copil de la Rebolledo-Jaramillo și colab. (2014). Fiecare culoare reprezintă un țesut diferit sau un proces de dezvoltare. Frunzele copacului reprezintă țesuturile epiteliale de sânge și obraz prelevate de la mamă și copilul ei. Liniile solide prezintă procese modelate de o cantitate fixă ​​de derivă genetică, iar liniile punctate prezintă procese în care deriva genetică se acumulează liniar cu vârsta. Componenta roșie, reprezentând oogeneza timpurie, modelează un blocaj de o singură generație, cu dublarea ulterioară a dimensiunii populației înapoi la o dimensiune mare. Descrierile parantetice în gri arată momentul evenimentelor de dezvoltare notabile.

Fiecare proces ontogenetic din filogenie este parametrizat de un parametru de derivare genetică și o rată de mutație. Rata mutației este unde este dimensiunea efectivă a populației de celule relevante și μ este rata de mutație per replicare, per bază. Deriva genetică poate fi modelată într-unul din cele trei moduri, și anume ca o cantitate fixă ​​de derivă genetică, ca o dimensiune explicită a blocajului sau ca o rată de acumulare a derivei genetice pe an.

Calculul probabilității

Dat copac ontogenetic cu k procesele ontogenetice, parametrii derivați genetici și ratele de mutație probabilitatea noastră este (1) unde reprezintă datele frecvenței heteroplasmelor. (Mai jos, indicele este lăsat pentru scurtă durată.) Să presupunem că frecvențele heteroplasmice au fost prelevate din F familii. Scrierea numărului de site-uri heteroplasmatice din familie eu, pentru datele privind frecvența heteroplasmelor la jlocusul heteroplasmatic (de) în familie eu, și pentru evenimentul respectiv j este heteroplasmatic în familie eu, probabilitatea noastră poate fi scrisă (2) unde este probabilitatea apariției variantelor heteroplasmatice în familie eu și este probabilitatea datelor heteroplasmice observate la jLocusul heteroplasmatic din familie eu, condiționată de heteroplasmie (adică, polimorfism) în cel puțin un țesut la acel locus. Presupunem că Poisson este distribuit cu rata unde G este dimensiunea genomului și este probabilitatea ca acest sit j este heteroplasmatic în familie eu. Observăm că pentru fiecare j și k adică probabilitatea heteroplasmiei depinde de familie (în mod specific, de vârsta indivizilor din familie) și nu de locusul particular.

Penalizăm partea din probabilitatea care implică numărul de variante heteroplasmatice cu parametrul α pentru a face inferența mai puțin sensibilă la detectarea experimentală a heteroplasmiei, care este o problemă netrivială, în special pentru heteroplasmii care segregează la frecvență joasă (Li și Stoneking 2012 Rebolledo-Jaramillo și colab. 2014). Fără o astfel de pedeapsă, probabilitatea este prea puternic influențată de numărul de heteroplasmii observate, o cantitate influențată atât de falsi pozitivi - cu o rată de până la heteroplasmele cu frecvență joasă în Rebolledo-Jaramillo și colab. (2014) - și de negative negative cauzate de praguri minime conservatoare ale frecvenței alelelor (1% în Rebolledo-Jaramillo și colab. 2014). Pe de altă parte, dacă numărul de heteroplasmii este complet absent din probabilitate, astfel încât toate informațiile despre derivă și mutație sunt luate numai din frecvențele heteroplasmice, distribuțiile posterioare ale ratelor de mutație sunt sensibile la frecvențele alelei exterioare care nu se potrivesc unui model. a derivei genetice și a mutației (rare), de asemenea. Ca un compromis, am stabilit valoarea acestei penalități de probabilitate la care, de fapt, se reduce artificial numărul total de site-uri luate în considerare în această componentă a probabilității, astfel încât, dacă, în realitate, sunt observate 500 de situri heteroplasmatice dintr-un total de site-uri, contribuția la probabilitate ar fi aceeași ca și dacă s-ar observa cinci situri heteroplasmatice într-un total de 1000 de situri.

Cu probabilitatea noastră (2), ignorăm implicit legătura dintre siturile heteroplasmatice dintr-o familie, chiar dacă, în realitate, lipsa recombinării înseamnă că site-urile sunt perfect legate. Justificăm această aproximare în două moduri: în primul rând, există de obicei puține variante heteroplasmatice care se regrupează într-o familie [înseamnă 2,6 în Rebolledo-Jaramillo și colab. (2014), 1,0 în Li și colab. (2016)] și, în al doilea rând, printre variantele heteroplasmatice care se regrupează într-o familie, cele mai multe segregează la frecvență scăzută, astfel încât modificările frecvenței unei heteroplasme nu afectează foarte mult frecvența alteia. Astfel, dinamica la mai multe situri heteroplasmatice ar trebui să semene cu cele ale unui model în care fiecare sit să separe cu adevărat independent. Această ipoteză este susținută de simulări ale genomurilor mitocondriale nerecombinante (vezi secțiunea Simulare de mai jos). Mai mult, presupunem că frecvențele heteroplasmelor sunt independente între familii.

Se determină că un sit este heteroplasmatic în conformitate cu etapele de filtrare descrise în Rebolledo-Jaramillo și colab. (2014), care includ filtre pentru calitatea cartografierii, calitatea bazei, frecvența minimă a alelei (), acoperirea (), complexitatea secvenței locale și contaminarea. Mai degrabă decât să calculăm probabilitățile pe baza frecvențelor de alele numite, modelăm eroarea de eșantionare binomială în numărul de citiri consensuale și alternative eșantionate dintr-o frecvență de alelă adevărată, necunoscută. Astfel reprezintă numărul de alele consens și alternative la jLocusul heteroplasmatic din familie eu. Condiționată de heteroplasmie (adică, polimorfism), probabilitatea citirii observate contează la locus j în familie eu este (3) unde este frecvența alelei adevărată, necunoscută, la locus j în familie eu. Suma se efectuează pe toate frecvențele alele posibile în țesuturile prelevate. Atât numeratorul, cât și numitorul pot fi calculate utilizând algoritmul de tăiere al lui Felsenstein (1981) - un algoritm de programare dinamică utilizat frecvent în calculele de probabilitate pentru arborii filogenetici. Detalii despre modul în care am calculat aceste cantități sunt date în apendicele A. Algoritmul de tăiere necesită calcularea distribuțiilor de tranziție a frecvenței alelelor pentru diferite valori ale parametrilor de derivare genetică și mutație. Așa cum este descris în Anexa B, am realizat acest lucru prin precomputarea numerică a distribuțiilor de tranziție a frecvenței alelelor sub modelul Wright-Fisher de generație discretă și interpolare liniară între valorile precomputate. Algoritmul de tăiere necesită, de asemenea, o distribuție a frecvențelor alelelor la rădăcina filogeniei, care, în aplicația noastră (vezi mai jos), reprezintă distribuția neobservabilă a frecvențelor alelelor heteroplasmice la mamă ca embrion. În urma lui Tătaru și colab. (2015), folosim o distribuție beta discretizată, simetrică, cu greutăți de probabilitate suplimentare, simetrice, la frecvențele 0 și 1. Cei doi parametri care specifică această distribuție sunt deduși în comun cu parametrii de derivare genetică și mutație.

Inferință

Luăm o abordare bayesiană a inferenței. Distribuțiile anterioare sunt Log-Uniform pentru parametrii genetici de derivare, măsurați în generații pe (de acum înainte „unități de derivare”). Pentru parametrii de deriva genetici specificați de o rată de acumulare a unităților de deriva pe an, limita inferioară (respectiv superioară) a limitelor de distribuție anterioare (uniforme) sunt împărțite la minimul (resp. Maxim) al vârstelor la care rata este multiplicat. Nu am permis ca efectele derivei genetice să scadă odată cu înaintarea în vârstă. Distribuțiile anterioare pentru dimensiunile blocajului sunt și, pentru parametrii ratei mutației, distribuția anterioară este Log-Uniform

Folosim o procedură MCMC de ansamblu afinar-invariant (Goodman și Weare 2010) pentru a preleva din distribuții posterioare, așa cum este implementat în pachetul Python emcee (Foreman-Mackey și colab. 2013). Evaluăm convergența prin inspecție vizuală a urmelor posterioare. Rularea a 500 de lanțuri în ansamblul MCMC pentru câte 20.000 de iterații fiecare, găsim o convergență bună după iterații și, astfel, aruncăm primele 5000 de iterații ale fiecărui lanț ca burn-in. Cu loci heteroplasmatici, o rulare durează 60-80 ore CPU, dar, datorită naturii paralele a ansamblului MCMC, calculele pot fi răspândite în mod eficient între procesoare, astfel încât, pe un nod de calcul cu 20 de nuclee, rezultatele sunt obținute în aproximativ 4 ore . Intervalele credibile (IC) raportate la 95% sunt intervale cu cea mai mare densitate posterioară.

Ca modalitate de evaluare a suportului relativ pentru diferite modele ontogenetice, estimăm factorii Bayes (adică, rapoartele integrale ale dovezilor posterioare) pentru modele ontogenetice alternative ale acumulării de derivă în linii celulare. Pentru modele și factorul Bayes este (4) unde este distribuția anterioară și este probabilitatea în cadrul modelului k. Aceste integrale de evidență posterioară sunt aproximate folosind emcee’s (Foreman-Mackey și colab. 2013) implementarea unei abordări care utilizează integrarea termodinamică (vezi Goggans și Chi 2004).

Simulare

Am efectuat două seturi de simulări pentru a testa procedura de inferență. Primele simulări au fost efectuate după modelul asumat de procedura noastră de inferență. Așa cum s-a descris mai sus, acest model presupune că fiecare locus se separă independent, tranzițiile de frecvență ale alelelor apar în conformitate cu modelul Wright-Fisher al derivei genetice și al mutației bialelelice, iar frecvențele heteroplasmei din rădăcina filogeniei ontogenetice sunt controlate de cei doi parametri ai unui discretizat. , distribuție beta simetrică cu greutate de probabilitate suplimentară la frecvențele zero și una. Aceste simulări au fost efectuate înainte în timp folosind un script Python personalizat.

Al doilea set de simulări a testat modul în care presupunerea noastră că segregarea loci afectează în mod independent inferența atunci când datele sunt simulate din genomi care nu sunt recombinați eșantionați din mai multe familii diferite. Aceste simulări au fost efectuate folosind o interfață personalizată pentru pachetul de simulare msprime (Kelleher și colab. 2016), care simulează variația genetică în cadrul modelului coalescent neutru standard cu mutație pe site-uri infinite. In these simulations, population sizes and branch lengths are equivalent to those under the forward-time simulations, but at the root of the ontogenetic phylogeny, we assume that ancestral lineages trace their ancestry back in time in a single panmictic population of constant size. Simulations were performed under conditions in which the distribution of the number of heteroplasmic variants per family roughly matched the distribution observed in the data.

We applied our inference procedure to a publicly available dataset, containing allele frequencies for 98 heteroplasmies sampled from 39 mother-offspring duos, originally published by Rebolledo-Jaramillo și colab. (2014). In this dataset, mitochondria from blood and cheek epithelial cells were sampled from both mother and offspring, resulting in a ontogenetic phylogeny with four leaves, each representing one of the four tissues sampled from a mother-offspring duo. Details of heteroplasmy discovery are described in Rebolledo-Jaramillo și colab. (2014).

To model the segregation of heteroplasmy frequencies during the ontogeny of the four tissues sampled from each duo, we used the ontogenetic phylogeny shown in Figure 1. This ontogenetic phylogeny models several life stages. The root of the phylogeny occurs at the divergence of the mother’s somatic and germline tissues when she is an embryo. On the branch leading to the somatic tissues in the mother, there is a brief period of early embryonic development before the blood and cheek epithelial cell lineages diverge at gastrulation as members of the ectodermal (cheek epithelial) and mesodermal (blood) germ layers. After diverging at gastrulation, each somatic tissue undergoes independent periods of genetic drift and mutation during later embryogenesis and early growth, and, finally, for each tissue there are independent rates of accumulation of genetic drift and mutation throughout adult life.

On the branch leading to the offspring tissues in the ontogenetic phylogeny in Figure 1, the first stage represented is the period of oogenesis prior to the birth of the mother, when the oogenic bottleneck is thought to occur. This is followed by the oocyte stage, during which we assume the mitochondria accumulate genetic drift and mutation at some rate linearly with the age of the mother before childbirth. At fertilization, this branch undergoes the same period of early somatic development experienced by the mother’s somatic tissues prior to gastrulation. Finally, the two somatic tissues of the offspring diverge at gastrulation and go through the same stages of development as the somatic tissues of the mother. For an overview of the events of human development, see, for example, Carlson (2014).

Effective oogenic bottleneck

Analyzing both simulated and real data, we find that there is limited power to infer the size of the oogenic bottleneck. This is to be expected, given that we also model the subsequent genetic drift of the later stages of oocyte development and in the early developing embryo each of these three ontogenetic processes occurs along the same branch of the ontogenetic phylogeny of the tissues considered here (Figure 1), which causes their respective contributions of genetic drift to be conflated with one another. We note that the genetic drift parameters of these ontogenetic processes are not truly unidentifiable: power to distinguish genetic drift during the early-oogenesis bottleneck from that of the later maternal germline is provided by the differing effects of genetic drift in mothers of different ages, and power to distinguish the contribution of drift in the early embryo is provided by the fact that this process occurs in both the mother and the offspring. Differences in effective population size (and thus scaled mutation rates) also provide theoretical power to distinguish these parameters, but, nevertheless, we find that these genetic drift parameters tend to become conflated with one another.

As a way of counteracting this conflation, we combine the genetic drift parameters of this branch in the ontogenetic phylogeny into an effective bottleneck size (EBS), summarizing the total genetic drift between mother and offspring. The effective bottleneck is comprised of the oogenic bottleneck în sine, the accumulation of genetic drift in the oocyte prior to ovulation, and the genetic drift in the embryo between fertilization and gastrulation. To combine genetic drift parameterized as a bottleneck with genetic drift parameterized in drift units, we used the approximate relationship where d is genetic drift in drift units, and is the bottleneck size. This approximation is justified in Appendix C. Using this relationship, our equation for the EBS has the form (5) where d is the summed genetic drift from the oogenic bottleneck în sine and pregastrulation embryogenesis, λ is the rate of genetic drift accumulation in the oocyte, and A is the age of the mother at childbirth. Because, in our model, genetic drift accumulates in the oocyte as the mother ages prior to ovulation, the size of the effective bottleneck decreases with age. We summarize this rate of decrease by linearizing (5) between ages 25 and 34, the first and third quartiles of maternal age at childbirth in the dataset from Rebolledo-Jaramillo și colab. (2014).

Data availability

Our inference procedure is released under a permissive license in a Python package called mope, available at https://github.com/ammodramus/mope or from the Python Package Index (PyPI, http://pypi.python.org/). As we describe above, our inference procedure requires precomputed transition distributions. These can be generated by the user or downloaded from https://github.com/ammodramus/mope. Our simulation scripts are also provided with the inference procedure.

Data from Rebolledo-Jaramillo și colab. (2014) are available from that paper’s supplemental material and from the NCBI Sequence Read Archive (www.ncbi.nlm.nih.gov/sra), accession SRP047378.


  • Boitard et al. (2016) Boitard, S., W. Rodríguez, F. Jay, S. Mona, and F. Austerlitz. 2016. Inferring population size history from large samples of genome-wide molecular data-an approximate Bayesian computation approach. PLoS genetics 12:e1005877.
  • Bollback et al. (2008) Bollback, J. P., T. L. York, and R. Nielsen. 2008. Estimation of 2Nes from temporal allele frequency data. Genetics 179:497–502.
  • Bouckaert (2010) Bouckaert, R. R. 2010. Densitree: making sense of sets of phylogenetic trees. Bioinformatics 26:1372–1373.
  • Bouckaert et al. (2019) Bouckaert, R., T. G. Vaughan, J. Barido-Sottani, S. Duchêne, M. Fourment, A. Gavryushkina, J. Heled, G. Jones, D. Kühnert, N. De Maio and others. 2019. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS computational biology 15, 4:e1006650.
  • Bryant et al. (2012) Bryant, D., R. Bouckaert, J. Felsenstein, N. A. Rosenberg, and A. RoyChoudhury. 2012. Inferring Species Trees Directly from Biallelic Genetic Markers: Bypassing Gene Trees in a Full Coalescent Analysis. Molecular Biology and Evolution 29:1917–1932.
  • Bryant and Hahn (2019) Bryant, D. and M. Hahn. 2019. The concatenation question. în Phylogenetics in the Genomics Era (F. Delsuc, N. Galtier, and C. Scornavacca, eds.). (in presa).
  • Chifman and Kubatko (2014) Chifman, J. and L. Kubatko. 2014. Quartet inference from SNP data under the coalescent model. Bioinformatics 30:3317–3324.
  • Drummond and Rambaut (2007) Drummond, A. J. and A. Rambaut. 2007. Beast: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology 7:214.
  • Edwards (2009) Edwards, S. V. 2009. Is a new and general theory of molecular systematics emerging? Evolution: International Journal of Organic Evolution 63:1–19.
  • Epstein and Mazzeo (2010) Epstein, C. L. and R. Mazzeo. 2010. Wright–Fisher diffusion in one dimension. SIAM Journal on Mathematical Analysis 42:568–608.
  • Etheridge (2011) Etheridge, A. 2011. Some mathematical models from population genetics : École d’éte de probabilités De Saint-flour XXXIX-2009. Springer.
  • Ethier and Norman (1977) Ethier, S. N. and M. F. Norman. 1977. Error estimate for the diffusion approximation of the Wright–Fisher model. Proceedings of the National Academy of Sciences 74:5096–5098.
  • Ewens (2004) Ewens, W. 2004. Mathematial population genetics. I. Theoretical introduction. Interdisciplinary Applied Mathematics 2 ed. Springer, New York.
  • Felsenstein (1981) Felsenstein, J. 1981. Evolutionary trees from dna sequences: a maximum likelihood approach. Journal of Molecular Evolution 17:368–376.
  • Felsenstein (1992) Felsenstein, J. 1992. Phylogenies from restriction sites: a maximum-likelihood approach. Evolution 46:159–173.
  • Fox and Parker (1968) Fox, L. and I. Parker. 1968. Chebyshev Polynomials in Numerical Analysis. Oxford Mathematical Handbooks Oxford University Press, London.
  • Georges et al. (2018) Georges, A., B. Gruber, G. B. Pauly, D. White, M. Adams, M. J. Young, A. Kilian, X. Zhang, H. B. Shaffer, and P. J. Unmack. 2018. Genomewide SNP markers breathe new life into phylogeography and species delimitation for the problematic short-necked turtles (Chelidae: Emydura) of eastern Australia. Molecular Ecology 27:5195–5213.
  • Griffiths and Spano‘ (2010) Griffiths, R. C. and D. Spano‘. 2010. Diffusion processes and the coalescent. Pages 358–379 în Probability and Mathematical Genetics: Papers in Honour of Sir John Kingman (N. H. Bingham and C. M. Goldie, eds.) London Mathematical Society Lecture Note Series. Cambridge University Press.
  • Gutenkunst et al. (2010) Gutenkunst, R., R. Hernandez, S. Williamson, and C. Bustamante. 2010. Diffusion approximations for demographic inference: DaDi. Nature Precedings 5:1–1.
  • Gutenkunst et al. (2009) Gutenkunst, R. N., R. D. Hernandez, S. H. Williamson, and C. D. Bustamante. 2009. Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLoS Genetics 5.
  • Heled et al. (2013) Heled, J., D. Bryant, and A. J. Drummond. 2013. Simulating gene trees under the multispecies coalescent and time-dependent migration. BMC evolutionary biology 13:44.
  • Hiscott et al. (2016) Hiscott, G., C. Fox, M. Parry, and D. Bryant. 2016. Efficient Recycled Algorithms for Quantitative Trait Models on Phylogenies. Genome biology and evolution 8:1338–1350.
  • Jenkins et al. (2017) Jenkins, P. A., D. Spano, et al. 2017. Exact simulation of the Wright–Fisher diffusion. The Annals of Applied Probability 27:1478–1509.
  • Kingman (1982) Kingman, J. 1982. The coalescent. Stochastic Processes and their Applications 13:235–248.
  • Lapierre et al. (2017) Lapierre, M., A. Lambert, and G. Achaz. 2017. Accuracy of demographic inferences from the site frequency spectrum: the case of the Yoruba population. Genetics 206:439–449.
  • Larget et al. (2010) Larget, B. R., S. K. Kotha, C. N. Dewey, and C. Ané. 2010. BUCKy: gene tree/species tree reconciliation with Bayesian concordance analysis. Bioinformatics 26:2910–2911.
  • Liu (2008) Liu, L. 2008. Best: Bayesian estimation of species trees under the coalescent model. Bioinformatics 24:2542–2543.
  • Liu et al. (2008) Liu, L., D. K. Pearl, R. T. Brumfield, and S. V. Edwards. 2008. Estimating species trees using multiple-allele DNA sequence data. Evolution: International Journal of Organic Evolution 62:2080–2091.
  • Liu et al. (2010) Liu, L., L. Yu, and S. V. Edwards. 2010. A maximum pseudo-likelihood approach for estimating species trees under the coalescent model. BMC Evolutionary Biology 10:302.
  • Lukic and Hey (2012) Lukic, S. and J. Hey. 2012. Demographic inference using spectral methods on SNP data, with an analysis of the human out-of-Africa expansion. Genetics 192:619–639.
  • Mason and Handscomb (2002) Mason, J. C. and D. C. Handscomb. 2002. Chebyshev polynomials. Chapman and Hall/CRC.
  • McKane and Waxman (2007) McKane, A. J. and D. Waxman. 2007. Singular solutions of the diffusion equation of population genetics. Journal of Theoretical Biology 247:849–858.
  • Nielsen (2000) Nielsen, R. 2000. Estimation of population parameters and recombination rates from single nucleotide polymorphisms. Genetics 154:931–942.
  • Øksendal (2003) Øksendal, B. 2003. Stochastic differential equations, an Introduction with Applications. Universitext Springer.
  • Racimo et al. (2016) Racimo, F., G. Renaud, and M. Slatkin. 2016. Joint estimation of contamination, error and demography for nuclear dna from ancient humans. PLoS genetics 12:e1005972.
  • Rambaut et al. (2018) Rambaut, A., A. J. Drummond, D. Xie, G. Baele, and M. A. Suchard. 2018. Posterior summarization in Bayesian phylogenetics using tracer 1.7. Systematic Biology 67:901–904.
  • Sirén et al. (2010) Sirén, J., P. Marttinen, and J. Corander. 2010. Reconstructing population histories from single nucleotide polymorphism data. Molecular Biology and Evolution 28:673–683.
  • Song and Steinrücken (2012) Song, Y. S. and M. Steinrücken. 2012. A simple method for finding explicit analytic transition densities of diffusion processes with general diploid selection. Genetics 190:1117–29.
  • Steinrücken et al. (2014) Steinrücken, M., A. Bhaskar, and Y. S. Song. 2014. A novel spectral method for inferring general diploid selection from time series genetic data. The Annals of Applied Statistics 8:2203–2222.
  • Tataru et al. (2017) Tataru, P., M. Simonsen, T. Bataillon, and A. Hobolth. 2017. Statistical inference in the Wright–Fisher model using allele frequency data. Systematic biology 66:e30–e46.
  • Trefethen (2013) Trefethen, L. N. 2013. Approximation theory and approximation practice. SIAM, Philadelphia.
  • Vachaspati and Warnow (2015) Vachaspati, P. and T. Warnow. 2015. ASTRID: accurate species trees from internode distances. BMC Genomics 16:S3.
  • Waldvogel (2006) Waldvogel, J. 2006. Fast construction of the Fejér and Clenshaw–Curtis quadrature rules. BIT Numerical Mathematics 46:195–202.
  • Williamson et al. (2005) Williamson, S. H., R. Hernandez, A. Fledel-Alon, L. Zhu, R. Nielsen, and C. D. Bustamante. 2005. Simultaneous inference of selection and population growth from patterns of variation in the human genome. Proceedings of the National Academy of Sciences 102:7882–7887.
  • Wright (1931) Wright, S. 1931. Evolution in mendelian populations. Genetics 16:97.
  • Yang (2015) Yang, Z. 2015. The BPP program for species tree estimation and species delimitation. Current Zoology 61:854–865.

A.1 The O ( k ) algorithm for solving diffusions

In the main text, we introduced a matrix Q which plays a key role in the numerical solution of diffusions. The matrix is chosen so that if

( β 1 ( 1 − x ) − β 2 x ) d d x g ( x , t ) + 1 2 x ( 1 − x ) d 2 d x 2 g ( x , t ) ≈ K ∑ k = 0 ( Q λ ( t ) ) k T ∗ k ( x ) . (23)

The bottleneck in the numerical method we use is the repeated solution of linear systems that look like

for difference complex values z and vectors v . Using a direct method, these take O ( K 2 ) time each as Q is lower triangular. However we can do better, using a trick described in Fox and Parker (1968) . Here we give a very high level description of the approach.

The key idea is to apply integration twice to the LHS of ( 23 ). Using integration by parts multiple times we have

∫ x 0 ∫ y 0 [ ( β 1 ( 1 − z ) − β 2 z ) d d z g ( z , t ) + 1 2 z ( 1 − z ) d 2 d z 2 g ( z , t ) ] d z d y
= − ∫ x 0 ∫ y 0 g ( z , t ) d z + ( β 1 ( 1 − x ) − β 2 x − 1 + 2 x ) ∫ x 0 g ( z , t ) d z
+ 1 2 x ( 1 − x ) g ( x , t ) + ( 1 2 − β 1 ) x g ( 0 , t ) . (24)

We introduce two new matrices X and Y . The matrix X is derived from properties of shifted Chebyshev polynomials, and is defined so that if g is expanded as in ( 22 ) then

The matrix Y comes from ( 24 ) and has the property that if g is expanded as in ( 22 ) then the RHS of ( 24 ) equals

The usefulness of this follows from that fact that, with the exception of two rows, all the non-zero entries in X and Y are on or near the digaonal: both matrices are almost banded. To solve ( Q − z I ) x = v we multiply both sides by X and solve the sparse system that results. Overall, this now takes O ( K ) time.

A.2 Further details on the simulation studies

Simulated from Priors θ 1 θ 2 θ 3 θ 4 θ 5 θ 6 θ r o o t height
Γ (2,200) Correct 0.91 0.91 0.94 0.93 0.94 0.97 0.95 1.00
Γ (2,200) Incorect 0.84 0.88 0.87 0.84 0.75 0.71 0.06 0.99
Γ (4,400) Correct 0.93 0.97 0.92 0.97 0.9 0.98 0.95 1.00
Γ (4,400) Incorect 0.91 0.95 0.93 0.94 0.78 0.68 0.32 0.99
Γ (80,8000) Correct 0.98 0.97 0.99 0.99 0.97 0.98 0.97 1.00
Γ (80,8000) Incorect 0.94 0.95 0.96 0.93 0.84 0.83 0.51 1.00
Table 3 : Frequency that parameters fall within the 95 % highest posterior density of the MCMC chain for a 100 simulated datasets varying the prior on the population sizes for the trees from which the data is simulated. Note that we display the frequencies above as percentages.

We conducted a simulation study to check whether the model posterior distributions recover model parameters. We simulate data from three different 4 species tree distributions. Parameter means are equal but we change the variance of the 4 species tree distribution. We then draw a 100 parameter sets from each 4 species tree distribution. Thereafter we simulate a 1000 SNPs for 32 individuals distributed over 4 species for each parameter set. For every simulated SNP dataset we ran two chains for 100,000 iterations, one with correct prior, i.e 4 species tree distribution used for simulation, and one with incorrect prior. See Table

4 for details on priors. We then count the number of parameters that were recovered, i.e fall within the 95% highest posterior density of the MCMC chain. MCMC chains took on average ∼ 200 seconds to run. We report these counts in Table 3 . We see that correct priors on the model has a high rate of parameter recovery. However, when incorrect priors are used recovery rates suffer for population size parameters close to the root. This does not come as a surprise due to information loss on branches as we go up the tree.