Informations

8.7 : Annexe - Algorithme d'élagage de Felsenstein - Biologie


L'algorithme d'élagage de Felsenstein (1973) est un exemple de programmation dynamique, un type d'algorithme qui a de nombreuses applications en biologie comparée. En programmation dynamique, nous décomposons un problème complexe en une série d'étapes plus simples qui ont une structure imbriquée. Cela nous permet de réutiliser les calculs de manière efficace et accélère le temps nécessaire pour faire des calculs.

La meilleure façon d'illustrer l'algorithme de Felsenstein est à travers un exemple, qui est présenté dans les panneaux ci-dessous. Nous essayons de calculer la probabilité d'un caractère à trois états sur un arbre phylogénétique qui comprend six espèces.

Figure 8.4A. Chaque pointe et nœud interne de l'arbre a trois cases, qui contiendront les probabilités pour les trois états de caractère à ce point de l'arbre. La première case représente un état 0, le deuxième état 1 et le troisième état 2. Image de l'auteur, peut être réutilisée sous licence CC-BY-4.0.

  1. La première étape de l'algorithme consiste à remplir les probabilités pour les pourboires. Dans ce cas, nous connaissons les états aux extrémités de l'arbre. Mathématiquement, nous affirmons que nous connaissons précisément les états des caractères aux extrémités ; la probabilité que cette espèce ait l'état que nous observons est 1, et tous les autres états ont une probabilité zéro :
  1. Ensuite, nous identifions un nœud où tous ses descendants immédiats sont des astuces. Il y aura toujours au moins un tel nœud ; souvent, il y en aura plus d'un, auquel cas nous en choisirons un arbitrairement. Pour cet exemple, nous choisirons le nœud qui est l'ancêtre commun le plus récent des espèces A et B, étiqueté nœud 1 sur la figure 8.2B.
  2. Nous utilisons ensuite l'équation 7.6 pour calculer la vraisemblance conditionnelle pour chaque état de caractère pour le sous-arbre qui inclut le nœud que nous avons choisi à l'étape 2 et ses descendants de pointe. Pour chaque état de caractère, la vraisemblance conditionnelle est la probabilité, compte tenu des données et du modèle, d'obtenir les états de caractère de pointe si vous commencez avec cet état de caractère à la racine. En d'autres termes, nous gardons une trace de la probabilité pour les parties de pointe de l'arbre, y compris nos données, si le nœud que nous considérons avait chacun des états de caractère possibles. Ce calcul est :

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

je et X sont les deux indices pour le k états des caractères, avec des sommes prises sur tous les états possibles aux extrémités des branches (X), et les termes calculés pour chaque état possible au nœud (je). Les deux éléments de l'équation sont les descendants gauche et droit du nœud d'intérêt. Les branches peuvent être attribuées arbitrairement à gauche ou à droite sans affecter le résultat final, et l'approche fonctionne également pour les polytomies (mais l'équation est légèrement différente). De plus, chacune de ces deux pièces a elle-même deux parties : la probabilité de début et de fin avec chaque état le long des deux branches considérées, et les probabilités conditionnelles actuelles qui entrent dans l'équation aux extrémités du sous-arbre (LL(X) et LR(X)). Les longueurs de branches sont désignées par tL et tR respectivement à gauche et à droite.

On peut penser à la vraisemblance "s'écoulant" le long des branches de l'arbre, et les vraisemblances conditionnelles pour les branches gauche et droite sont combinées via la multiplication à chaque nœud, générant la vraisemblance conditionnelle pour le nœud parent pour chaque état de caractère (LP(je)).

Considérons le sous-arbre menant aux espèces A et B dans l'exemple donné. Les deux états de caractère de pointe sont 0 (pour l'espèce A) et 1 (pour l'espèce B). Nous pouvons calculer la vraisemblance conditionnelle pour l'état de caractère 0 au nœud 1 comme :

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

Ensuite, nous pouvons calculer les termes de probabilité à partir de la matrice de probabilité P. Dans ce cas tL = tR = 1,0, donc pour les branches gauche et droite :

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

Pour que:

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

Notez maintenant que, puisque l'état du caractère gauche est connu pour être zéro, LL(0)=1 et LL(1)=LL(2)=0. De même, le bon état est un, donc LR(1)=1 et LR(0)=LR(2)=0.

Nous pouvons maintenant remplir les deux parties de l'équation 8.2 :

$$ sumlimits_{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} $$

et:

$$ sumlimits_{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} $ $

Donc:

LP(0)=0.37 ⋅ 0.32 = 0.12.

Cela signifie que sous le modèle, si l'état au nœud 1 était 0, nous aurions une probabilité de 0,12 pour cette petite section de l'arbre. Nous pouvons utiliser une approche similaire pour trouver que :

(éq. 8.13)

LP(1)=0.32 ⋅ 0.37 = 0.12.

LP(2)=0.32 ⋅ 0.32 = 0.10.

Maintenant, nous avons la probabilité pour les trois états ancestraux possibles. Ces numéros peuvent être entrés dans les cases appropriées :

  1. Nous répétons ensuite le calcul ci-dessus pour chaque nœud de l'arbre. Pour les nœuds 3 à 5, tous les LL(X) et LR(X) les termes sont nuls ; leurs valeurs peuvent être lues dans les cases de l'arbre. Le résultat de tous ces calculs :
  1. Nous pouvons maintenant calculer la vraisemblance sur l'ensemble de l'arbre en utilisant les vraisemblances conditionnelles pour les trois états à la racine de l'arbre.

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

??X est la probabilité a priori de cet état de caractère à la racine de l'arbre. Pour cet exemple, nous considérerons que ces probabilités a priori sont uniformes, égales pour chaque état (??X = 1/k = 1/3). La vraisemblance pour notre exemple est donc :

[L = 1/3 0,00150 + 1/3 0,00151 + 1/3 ⋅ 0,00150 = 0,00150 label{8.15}]

Notez que si vous essayez cet exemple dans un autre logiciel, comme GEIGER ou PAUP*, le logiciel calculera une probabilité ln de -6,5, qui est exactement le logarithme naturel de la valeur calculée ici.


Un problème de reconstruction pour une classe de réseaux phylogénétiques avec transferts de gènes latéraux

Les transferts de gènes latéraux ou horizontaux sont un type d'événements évolutifs asymétriques où le matériel génétique est transféré d'une espèce à une autre. Dans cet article, nous considérons Réseaux LGT, un modèle général de réseaux phylogénétiques avec des transferts de gènes latéraux qui consistent, grosso modo, en un principal arbre enraciné avec ses feuilles étiquetées sur un ensemble de taxons, et un ensemble de secondaire des arcs entre les nœuds de cet arbre représentant des transferts de gènes latéraux. Un réseau LGT donne naissance de manière naturelle à un sous-arbre phylogénétique principal et un ensemble de sous-arbres phylogénétiques secondaires, qui, grosso modo, représentent respectivement la ligne principale d'évolution de la plupart des gènes et les lignes secondaires d'évolution par transferts latéraux de gènes.

Résultats

Nous introduisons un ensemble de conditions simples sur un réseau LGT qui garantissent que ses sous-arbres phylogénétiques principal et secondaire sont deux à deux différents et que ces sous-arbres déterminent, à isomorphisme près, le réseau LGT. Nous donnons ensuite un algorithme qui, étant donné un ensemble d'arbres phylogénétiques différents deux à deux (T_0,T_1,ldots ,T_k) sur le même ensemble de taxons, sort, lorsqu'il existe, le réseau LGT qui satisfait ces conditions et tel que son arbre phylogénétique principal est (T_0) et ses arbres phylogénétiques secondaires sont (T_1,ldots ,T_k) .


Résumé

Les modèles de coévolution moléculaire peuvent révéler des contraintes structurelles et fonctionnelles au sein ou entre les molécules organiques. Ces modèles sont mieux compris lorsque l'on considère le processus évolutif sous-jacent, ce qui nous permet de démêler le signal de l'évolution dépendante des sites (coévolution) des effets de l'ascendance partagée des gènes. Inversement, le fait de ne pas tenir compte de l'évolution dépendante des sites lors de l'étude de l'histoire des gènes a un impact négatif sur la précision des arbres phylogénétiques inférés. Bien que la coévolution moléculaire et l'histoire phylogénétique soient interdépendantes, les analyses des deux processus sont menées séparément, un choix dicté par la commodité des calculs, mais au détriment de la précision. Nous présentons une méthode bayésienne et un logiciel associé pour déduire combien et quels sites d'un alignement évoluent selon un processus évolutif indépendant ou dépendant par paires, et pour estimer simultanément les relations phylogénétiques entre les séquences. Nous validons notre méthode sur des ensembles de données synthétiques et défions nos prédictions de coévolution sur la molécule d'ARNr 16S en les comparant avec sa structure moléculaire connue. Enfin, nous évaluons l'exactitude des arbres phylogénétiques déduits sous l'hypothèse d'indépendance entre les sites en utilisant des ensembles de données synthétiques, la molécule d'ARNr 16S et 10 alignements supplémentaires de gènes codant pour des protéines d'eucaryotes. Nos résultats démontrent qu'inférer des arbres phylogénétiques tout en tenant compte de l'évolution du site dépendant a un impact significatif sur les estimations de la phylogénie et du processus évolutif.

La coévolution moléculaire est le processus évolutif par lequel les interactions entre des sites distants d'une molécule, ou des sites de molécules différentes, sont maintenues de manière à préserver des propriétés fonctionnelles ou structurelles avantageuses. Par exemple, des fragments coévoluant au sein de séquences protéiques sont impliqués dans les contraintes de repliement et informent sur les intermédiaires de repliement, l'assemblage de peptides ou des mutations clés avec des rôles connus dans les maladies génétiques (1, 2). La disponibilité toujours croissante de séquences moléculaires (nucléotides et acides aminés) nous fournit une quantité sans précédent de données qui ont un fort potentiel pour révéler des gènes et des régions de gènes évoluant selon un processus contraint (3, 4).

Il existe plusieurs méthodes pour déduire la coévolution à partir des données de séquence uniquement (basées sur des modèles d'appariement entre les sites, comme examiné dans les références 5 et 6). Cependant, ces méthodes n'exploitent pas un élément clé dans la modélisation des processus évolutifs sous-jacents : l'arbre phylogénétique décrivant les relations entre les séquences moléculaires. L'intégration du signal phylogénétique dans l'analyse de la coévolution est cruciale car elle permet de distinguer les motifs véritablement coévolutifs des motifs similaires induits par l'histoire commune des séquences (5, 7). À cette fin, plusieurs méthodes ont été développées pour déduire la coévolution tout en tenant compte des relations phylogénétiques (7), mais seules quelques-unes d'entre elles modélisent explicitement le processus de coévolution le long d'un arbre phylogénétique donné (8 ⇓ ⇓ –11).

Toutes les méthodes sensibles à la phylogénie pour détecter la coévolution reposent sur l'hypothèse que les relations phylogénétiques entre les séquences sont connues et peuvent être traitées comme des « données observées ». Typiquement, les arbres phylogénétiques sont eux-mêmes inférés à partir de données moléculaires (12), mais leur inférence est basée sur une hypothèse fondamentale selon laquelle chaque site évolue indépendamment de tous les autres (13). Cette hypothèse, qui est évidemment violée en présence de coévolution, présente des avantages en termes de traçabilité informatique, car la probabilité d'un alignement étant donné un arbre phylogénétique est le produit de la probabilité individuelle de chaque site. Il a été démontré que cette simplification du mécanisme évolutif en présence de sites non indépendants diminue la précision des arbres phylogénétiques inférés (14, 15). Cependant, les ensembles de données avec de fortes contraintes fonctionnelles ou structurelles sont souvent analysés dans des cadres phylogénétiques qui supposent l'indépendance entre les sites. Par exemple, la petite sous-unité ribosomique (16S) est fréquemment utilisée pour estimer les premières relations évolutives entre les principales lignées de l'arbre de vie (16, 17), négligeant ses nombreuses contraintes structurelles et les preuves de coévolution (18).

La présence de motifs coévolutifs dans de nombreuses séquences de nucléotides et d'acides aminés s'étend bien au-delà du gène 16S et est étayée par un grand nombre de preuves (9, 19). Ignorer les interdépendances entre l'histoire phylogénétique des séquences et les processus contraints régissant l'évolution des nucléotides ou des acides aminés peut gravement entraver notre capacité à déduire des arbres phylogénétiques corrects (14, 15) et à détecter avec précision la coévolution (5, 7). Cependant, l'inférence de ces processus interdépendants est toujours menée séparément pour des raisons mathématiques, au détriment de l'hypothèse d'une indépendance entre les sites qui a été initialement décrite comme « pas nécessairement biologiquement valide » par Felsenstein en 1983, à l'aube de la phylogénétique moléculaire basée sur la vraisemblance ( 13). Pour résoudre ce problème, nous présentons un cadre bayésien pour analyser un alignement de nucléotides et estimer conjointement (je) le nombre de paires de sites qui coévoluent et leur position dans la séquence, se différenciant ainsi d'un modèle de fond d'évolution indépendante, (ii) les paramètres des modèles de substitution indépendants et dépendants, et (iii) l'arbre phylogénétique décrivant les relations entre les séquences. Notre méthode s'appelle CoevRJ et le logiciel qui l'implémente est disponible sur https://bitbucket.org/XavMeyer/coevrj.

Nous évaluons les performances de CoevRJ dans la reconstruction d'arbres phylogénétiques et l'inférence de coévolution sur la base d'une vaste gamme d'ensembles de données simulés, d'un alignement de l'ARNr 16S hautement coévolutif et de 10 ensembles de données empiriques eucaryotes de gènes codant des protéines. Nous montrons que CoevRJ fournit une identification précise des sites coévoluants, validée par des données simulées et empiriques. Nous évaluons les effets de la coévolution sur les estimations phylogénétiques en comparant nos résultats avec ceux obtenus sous l'hypothèse d'une évolution indépendante et démontrons l'importance de tenir compte de la dépendance entre les sites lors de l'inférence des arbres phylogénétiques sur les ensembles de données soumis à la coévolution.


Méthodes basées sur la distance

Après avoir lu l'alignement, nous pouvons construire un premier arbre avec des méthodes basées sur la distance. La fonction dist.dna du package ape calcule les distances pour de nombreux modèles de substitution d'ADN. Pour utiliser la fonction dist.dna, nous devons transformer les données en classe DNAbin. Pour les acides aminés, la fonction dist.ml propose des modèles de substitution courants (par exemple "WAG", "JTT", "LG", "Dayhoff", "cpREV", "mtmam", "mtArt", "MtZoa" ou "mtREV24" ).

Après avoir construit une matrice de distance, nous reconstruisons un arbre enraciné avec UPGMA et alternativement un arbre non enraciné en utilisant Neighbor Joining (Saitou et Nei 1987) (Studier et Keppler 1988) . Plus de méthodes de distance comme fastme sont disponibles dans le singe emballer.

On peut tracer les arbres treeUPGMA et treeNJ avec les commandes :

Les méthodes basées sur la distance sont très rapides et nous utiliserons les arbres UPGMA et NJ comme arbres de départ pour les analyses de parcimonie maximale et de vraisemblance maximale.


Conclusion

Le séquençage d'un grand nombre de génomes entiers a permis d'étudier les schémas de gain et de perte de gènes à une échelle énorme. Même si les méthodes pour déduire les gains et les pertes existent depuis près de 30 ans, ce n'est qu'avec l'analyse de génomes entiers que l'effet de petits biais devient clair. L'effet net du biais dans les méthodes de réconciliation des arbres démontré ici est que le nombre de pertes de gènes déduites ne doit pas être pris pour argent comptant, et le nombre de duplications déduites doit être analysé avec soin.

Comment surmonter le biais des méthodes de réconciliation actuelles ? Les algorithmes qui fournissent une estimation du support statistique pour chaque gain ou perte inféré (par exemple, [12, 30]) ou qui prennent en compte la longueur des branches d'arbres d'espèces peuvent tous deux offrir des améliorations aux méthodes actuelles. Cette dernière possibilité offre un moyen d'améliorer les inférences car les branches d'arbres d'espèces courtes sont celles qui sont les plus susceptibles de conduire à des arbres de gènes inférés à tort, comme on le voit dans le cas d'un tri incomplet des lignées. Il est également important de noter que la plupart des biais décrits ici se produisent lorsqu'il y a un nombre égal de gènes parmi les taxons - s'il y a un nombre inégal de gènes, alors des duplications et des pertes peuvent être déduites des informations de présence/absence. Dans ce cas, le problème se réduit simplement à un problème concernant l'évolution du nombre de copies, ce qui est exactement l'approche de la méthode de vraisemblance mentionnée ci-dessus. Mais cette méthode n'est pas non plus sans biais [3].

Enfin, les biais dans les méthodes de réconciliation des arbres jettent le doute sur les travaux antérieurs sur l'évolution des génomes des vertébrés. Cependant, les résultats présentés ici ne peuvent pas réfuter les résultats précédents, car un excès d'anciens doublons et des pertes de gènes récentes ont été déduits même en utilisant des seuils bootstrap de 100 %. La découverte du véritable modèle d'évolution du génome des vertébrés nécessitera des résultats de simulation, de meilleurs arbres génétiques, plus de données ou une combinaison des trois.


Conclusion

Nous introduisons le MowgliNNI méthode heuristique s'appuyant sur des réarrangements NNI des parties incertaines des arbres génétiques pour résoudre un problème d'optimisation de parcimonie pour les rapprochements tenant compte des duplications (), pertes () et transferts (). Nous montrons des preuves expérimentales que les rapprochements calculés selon le critère de parcimonie peuvent corriger efficacement les parties erronées des arbres de gènes déduits des données de séquence. Sur des données simulées, MowgliNNI propose souvent une nouvelle topologie d'arbre de gènes qui est plus proche de la bonne et qui conduit également à de meilleures prédictions. De plus, le nombre d'événements et le nombre de rapprochements les plus parcimonieux prédits par MowgliNNI sont nettement inférieurs à ceux obtenus sans remettre en cause la topologie de l'arbre génétique. Ceci est confirmé sur des données réelles. Un point critique pour les méthodes de parcimonie est le choix des coûts respectifs pour les événements évolutifs considérés. Nous montrons ici que MowgliNNILes performances de s ne sont que légèrement modifiées lors de la modification des coûts attribués aux événements individuels (, et), c'est-à-dire que la méthode est robuste aux erreurs de spécification des coûts.


Estimation du soutien

Étant donné T X et T oui dans une classe d'équivalence, une jointure T X oui sur paire commandée (T X,T oui) est un FST s'il est soutenu par au moins une fraction F des arbres d'entrée (c'est-à-dire si au moins une fraction F des arbres d'entrée ont T X oui comme sous-arbre). Notez que de tels T X oui est pris en charge uniquement par les arbres qui prennent en charge à la fois T X et T oui également. Motivé par cela, pour chaque FST T X nous maintenons un liste de soutien, noté par T X.supList, qui contient toutes les arborescences de la collection d'entrée qui prennent en charge T X. Pour estimer si la jointure sur (T X,T oui) aboutit à un FST, nous appliquons le théorème 1 uniquement sur les arbres dans T X.supList∩T oui.supListe. Nous stockons la liste de support sous forme de bitmap [36] pour une utilisation efficace de la mémoire et un calcul rapide de l'intersection des listes de support.


8.7 : Annexe - Algorithme d'élagage de Felsenstein - Biologie

Partie I : ANALYSE DE SÉQUENCE.

1. Alignement global par paires de séquences.

1.1 Alignement et évolution.

1.3 Un barème de notation pour le modèle.

1.4 Trouver les alignements les plus performants avec la programmation dynamique.

1.4.1 Déterminer Hje,j.

1.4.3 Trouver les alignements qui donnent le score le plus élevé.

1.6 Écarts de notation : pénalités d'écart.

1.7 Programmation dynamique pour la pénalité d'écart général.

1.8 Programmation dynamique pour la pénalité Affine Gap.

1.9 Score d'alignement et distance de séquence.

2 Alignement local par paire et recherche dans la base de données.

2.1 L'opération de base : comparer deux séquences.

2.3.2 Trouver les meilleurs alignements locaux.

2.3.4 Matrices de notation et pénalités d'écart.

2.4.2 Prétraiter la requête : faire la liste de mots.

2.4.3 Parcourir les séquences de la base de données.

3. Analyse statistique.

3.1 Test d'hypothèse pour l'homologie de séquence.

3.1.1 Génération aléatoire de séquences.

3.1.2 Utilisation de Z valeurs pour estimer la signification statistique.


Méthodes

Phylogénies ontogénétiques

Nous modélisons les mitochondries dans les tissus prélevés sur un ou plusieurs individus apparentés en tant que groupe de populations apparentées par une phylogénie ontogénétique. Le long de chaque branche de la phylogénie ontogénétique, les fréquences d'hétéroplasmie dans certains tissus ancestraux changent en raison de l'action de la dérive génétique et de la mutation. Nous supposons que la forme de la phylogénie ontogénétique est donnée.

Notre modèle de phylogénie ontogénétique diffère de plusieurs manières importantes du cadre typique de vraisemblance population-phylogénétique. Dans le modèle typique de génétique des populations, chaque branche est considérée comme une période indépendante de l'histoire de l'évolution et est donc sous le contrôle d'un paramètre indépendant. Contrairement à cela, nous permettons à un seul paramètre de déterminer la dynamique sur plusieurs parties de la phylogénie, car un seul processus de développement peut agir chez plusieurs individus apparentés, et ce processus de développement peut être supposé agir de la même manière chez chaque individu. De plus, alors qu'il est généralement supposé que chaque locus a été transmis par une seule phylogénie et a donc été soumis aux mêmes forces génétiques de la population, nous permettons aux effets de la dérive génétique et de la mutation de dépendre de l'âge des individus échantillonnés. En particulier, pour certains processus ontogénétiques, nous modélisons le taux d'accumulation de dérive génétique et de mutation avec l'âge. Ceci est motivé par des observations antérieures selon lesquelles les variants hétéroplasmiques se séparent et s'accumulent avec le temps dans les tissus somatiques (Sondheimer et al. 2011 Rebolledo-Jaramillo et al. 2014 Li et al. 2015) et au sein de la lignée germinale (Rebolledo-Jaramillo et al. 2014 Li et al. 2016 Wachsmuth et al. 2016). Enfin, dans le modèle phylogénétique de population typique, chaque branche de la phylogénie représente une seule période de l'histoire de l'évolution et est modélisée par un seul paramètre. Étant donné que plusieurs processus ontogénétiques d'intérêt peuvent se produire le long d'une seule branche d'une phylogénie ontogénétique, nous permettons aux branches de la phylogénie ontogénétique d'être divisées en plusieurs processus ontogénétiques distincts, contrôlés par des paramètres indépendants. La figure 1 montre ces caractéristiques avec une phylogénie ontogénétique représentant les relations entre deux tissus échantillonnés à la fois chez une mère et sa progéniture.

Phylogénie ontogénétique des tissus prélevés dans les duos mère-enfant de Rebolledo-Jaramillo et al. (2014). Chaque couleur représente un tissu ou un processus de développement différent. Les feuilles de l'arbre représentent le sang et les tissus épithéliaux des joues prélevés sur la mère et son enfant. Les lignes pleines montrent les processus modélisés par une quantité fixe de dérive génétique et les lignes en pointillés montrent les processus dans lesquels la dérive génétique s'accumule linéairement avec l'âge. La composante rouge, représentant l'ovogenèse précoce, modélise un goulot d'étranglement à génération unique avec un doublement ultérieur de la taille de la population jusqu'à une grande taille. Les descriptions entre parenthèses en gris montrent le moment des événements développementaux notables.

Chaque processus ontogénétique de la phylogénie est paramétré par un paramètre de dérive génétique et un taux de mutation. Le taux de mutation est où est la taille effective de la population cellulaire pertinente et ?? est le taux de mutation par réplication et par base. La dérive génétique peut être modélisée de l'une des trois manières suivantes, à savoir comme une quantité fixe de dérive génétique, comme une taille de goulot d'étranglement explicite ou comme un taux d'accumulation de dérive génétique par an.

Calcul de probabilité

Arbre ontogénétique donné avec k processus ontogénétiques, paramètres de dérive génétique et taux de mutation, notre probabilité est (1) où représente les données de fréquence d'hétéroplasmie. (Ci-dessous, l'indice est laissé de côté par souci de concision.) Supposons que les fréquences d'hétéroplasmie ont été échantillonnées à partir de F des familles. Ecrire pour le nombre de sites hétéroplasmiques dans la famille je, pour les données de fréquence d'hétéroplasmie au je locus hétéroplasmique (de ) dans la famille je, et pour l'événement ce site j est hétéroplasmique dans la famille je, notre vraisemblance peut s'écrire (2) où est la probabilité que des variants hétéroplasmiques surviennent dans la famille je et est la probabilité des données d'hétéroplasmie observées au je locus hétéroplasmique dans la famille je, conditionnelle à l'hétéroplasmie (c'est à dire., polymorphisme) dans au moins un tissu à ce locus. Nous supposons que Poisson est distribué avec un taux où g est la taille du génome, et est la probabilité que le site j est hétéroplasmique dans la famille je. Notons que pour chaque j et k c'est-à-dire que la probabilité d'hétéroplasmie dépend de la famille (en particulier de l'âge des individus de la famille) et non du locus particulier.

On pénalise la part de vraisemblance impliquant le nombre de variants hétéroplasmiques avec le paramètre ?? afin de rendre l'inférence moins sensible à la détection expérimentale de l'hétéroplasmie, qui est un problème non trivial, en particulier pour les hétéroplasmies ségrégeant à basse fréquence (Li et Stoneking 2012 Rebolledo-Jaramillo et al. 2014). Sans une telle pénalité, la vraisemblance est trop fortement influencée par le nombre d'hétéroplasmies observés, une quantité influencée à la fois par les faux positifs - à un taux allant jusqu'à pour les hétéroplasmies de basse fréquence à Rebolledo-Jaramillo et al. (2014)—et par les faux négatifs causés par les seuils de fréquences alléliques minimales conservatrices (1% à Rebolledo-Jaramillo et al. 2014). D'un autre côté, si le nombre d'hétéroplasmies est complètement absent de la probabilité, de sorte que toutes les informations sur la dérive et la mutation ne proviennent que des fréquences d'hétéroplasmie, les distributions postérieures des taux de mutation sont sensibles aux fréquences alléliques aberrantes qui ne correspondent pas à un modèle. de dérive génétique et de mutations (peu fréquentes). A titre de compromis, nous fixons la valeur de cette pénalité de vraisemblance à laquelle, en effet, réduit artificiellement le nombre total de sites considérés dans cette composante de la vraisemblance, de telle sorte que si, en réalité, 500 sites hétéroplasmiques sont observés sur un total de sites, la contribution à la probabilité serait la même que si cinq sites hétéroplasmiques étaient observés sur un total de 1000 sites.

Avec notre vraisemblance (2), nous ignorons implicitement la liaison entre les sites hétéroplasmiques au sein d'une famille, même si, en réalité, l'absence de recombinaison signifie que les sites sont parfaitement liés. Nous justifions cette approximation de deux manières : premièrement, il y a généralement peu de variants hétéroplasmiques coségrégeant dans une famille [moyenne 2,6 chez Rebolledo-Jaramillo et al. (2014), 1,0 dans Li et al. (2016)], et, deuxièmement, parmi les variants hétéroplasmiques coségrégeant dans une famille, la plupart se séparent à basse fréquence, de sorte que les changements de fréquence d'une hétéroplasmie n'affectent pas beaucoup la fréquence d'une autre. Ainsi, la dynamique de plusieurs sites hétéroplasmiques devrait ressembler étroitement à celle d'un modèle dans lequel chaque site se sépare vraiment indépendamment. Cette hypothèse est étayée par des simulations de génomes mitochondriaux non recombinants (voir section Simulation au dessous de). Nous supposons en outre que les fréquences d'hétéroplasmie sont indépendantes entre les familles.

Un site est déterminé comme étant hétéroplasmique selon les étapes de filtrage décrites dans Rebolledo-Jaramillo et al. (2014), qui incluent des filtres pour la qualité de la cartographie, la qualité de base, la fréquence allélique minimale ( ), la couverture ( ), la complexité de la séquence locale et la contamination. Plutôt que de calculer les probabilités en fonction des fréquences alléliques appelées, nous modélisons l'erreur d'échantillonnage binomiale dans le nombre de lectures consensuelles et alternatives échantillonnées à partir d'une fréquence allélique vraie et inconnue. Représente ainsi le nombre d'allèles consensus et alternatifs au je locus hétéroplasmique dans la famille je. Conditionnel à l'hétéroplasmie (c'est à dire., polymorphisme), la probabilité des comptes de lecture observés au locus j dans la famille je est (3) où est la fréquence allélique vraie, inconnue, au locus j dans la famille je. La somme est effectuée sur toutes les fréquences alléliques possibles dans les tissus échantillonnés. Le numérateur et le dénominateur peuvent être calculés à l'aide de l'algorithme d'élagage de Felsenstein (1981), un algorithme de programmation dynamique fréquemment utilisé dans les calculs de vraisemblance pour les arbres phylogénétiques. Des détails sur la façon dont nous avons calculé ces quantités sont donnés à l'annexe A. L'algorithme d'élagage nécessite le calcul des distributions de transition de fréquence allélique pour différentes valeurs de paramètres de dérive génétique et de mutation. Comme décrit dans l'annexe B, nous y sommes parvenus en précalculant numériquement les distributions de transition de fréquence allélique sous le modèle de Wright-Fisher à génération discrète, et en interpolant linéairement entre les valeurs précalculées. L'algorithme d'élagage nécessite également une distribution des fréquences alléliques à la racine de la phylogénie, qui, dans notre application (voir ci-dessous), représente la distribution inobservable des fréquences alléliques de l'hétéroplasmie chez la mère en tant qu'embryon. Suivre Tataru et al. (2015), nous utilisons une distribution bêta symétrique discrétisée avec des poids de probabilité supplémentaires et symétriques aux fréquences 0 et 1. Les deux paramètres spécifiant cette distribution sont inférés conjointement avec les paramètres de dérive génétique et de mutation.

Inférence

Nous adoptons une approche bayésienne de l'inférence. Les distributions antérieures sont Log-Uniformes pour les paramètres de dérive génétique, mesurés en générations par (ci-après « unités de dérive »). Pour les paramètres de dérive génétique spécifiés par un taux d'accumulation d'unités de dérive par an, la limite inférieure (resp. supérieure) des limites de distribution a priori (uniformes) est divisée par le minimum (resp. maximum) des âges auxquels le taux est multiplié. Nous n'avons pas laissé les effets de la dérive génétique diminuer avec l'âge. Les distributions antérieures sur les tailles des goulots d'étranglement sont et, pour les paramètres de taux de mutation, la distribution antérieure est Log-Uniforme

Nous utilisons une procédure MCMC d'ensemble affine-invariant (Goodman et Weare 2010) pour échantillonner à partir de distributions postérieures, comme implémenté dans le package Python emcee (Foreman-Mackey et al. 2013). Nous évaluons la convergence par inspection visuelle des traces postérieures. En exécutant 500 chaînes dans l'ensemble MCMC pour 20 000 itérations chacune, nous trouvons une bonne convergence après les itérations et, par conséquent, nous rejetons les 5000 premières itérations de chaque chaîne en tant que burn-in. Avec les loci hétéroplasmiques, une exécution prend 60 à 80 heures CPU, mais, en raison de la nature parallèle de l'ensemble MCMC, les calculs peuvent être efficacement répartis sur les CPU, de sorte que, sur un nœud de calcul à 20 cœurs, les résultats sont obtenus en ∼ 4 heures . Les intervalles de crédibilité (IC) à 95 % signalés sont des intervalles de densité postérieure la plus élevée.

Afin d'évaluer le support relatif de différents modèles ontogénétiques, nous estimons les facteurs de Bayes (c'est à dire., ratios d'intégrales de preuves postérieures) pour des modèles ontogénétiques alternatifs de l'accumulation de dérive au sein des lignées cellulaires. Pour les modèles et le facteur de Bayes est (4) où est la distribution a priori et est la vraisemblance sous le modèle k. Ces intégrales de preuves postérieures sont approximées à l'aide de l'animateur (Foreman-Mackey et al. 2013) mise en œuvre d'une approche utilisant l'intégration thermodynamique (voir Goggans et Chi 2004).

Simulation

Nous avons effectué deux séries de simulations pour tester notre procédure d'inférence. Les premières simulations ont été réalisées sous le modèle assumé par notre procédure d'inférence. Comme décrit ci-dessus, ce modèle suppose que chaque locus ségrége indépendamment, les transitions de fréquence allélique se produisent selon le modèle de Wright-Fisher de dérive génétique et de mutation biallélique, et les fréquences d'hétéroplasmie à la racine de la phylogénie ontogénétique sont contrôlées par les deux paramètres d'un modèle discrétisé. , distribution bêta symétrique avec un poids de probabilité supplémentaire aux fréquences zéro et un. Ces simulations ont été effectuées dans le temps à l'aide d'un script Python personnalisé.

La deuxième série de simulations a testé comment notre hypothèse selon laquelle les loci se séparent affecte indépendamment l'inférence lorsque les données sont simulées à partir de génomes non recombinants échantillonnés dans de nombreuses familles différentes. Ces simulations ont été effectuées à l'aide d'une interface personnalisée vers le package de simulation msprime (Kelleher et al. 2016), qui simule la variation génétique sous le modèle standard de coalescence neutre avec mutation en sites infinis. 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 et al. (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 et al. (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 en soi, 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 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 is the summed genetic drift from the oogenic bottleneck en soi and pregastrulation embryogenesis, ?? is the rate of genetic drift accumulation in the oocyte, and une 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 et al. (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 et al. (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. dans Phylogenetics in the Genomics Era (F. Delsuc, N. Galtier, and C. Scornavacca, eds.). (dans la presse).
  • 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 dans 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. La presse de l'Universite de Cambridge.
  • 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) Incorrect 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) Incorrect 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) Incorrect 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.