Module 5 · Section 5.4

Maximum de vraisemblance

60 min · niveau avance · dernière révision : 2026-11-01

Objectifs d'apprentissage

Écrire la vraisemblance d'un alignement sous un modèle d'évolution et la calculer par élagage
Distinguer paramètres de modèle et longueurs de branches, puis interpréter une log-vraisemblance
Expliquer la recherche d'arbres par NNI et SPR et choisir un modèle avec ModelFinder et le BIC
Interpréter des supports UFBoot et bootstrap sans les confondre avec une probabilité de vérité
Exécuter et documenter une analyse IQ-TREE ou RAxML-NG sur un alignement donné

Prérequis : 5.3 Neighbor-Joining et bootstrap

IQ-TREE 2.3.6RAxML-NG 1.2.0ModelFinder inclus dans IQ-TREE

Prérequis

Cette section prolonge la construction d'arbres par Neighbor-Joining et l'évaluation de leur robustesse vus en 5.3. Elle suppose connus les modèles d'évolution, les distances corrigées et la saturation présentés en 5.2, ainsi que la préparation d'un alignement multiple en 5.1. Les matrices de substitution et la distinction entre score et probabilité sont introduites en 2.5. Nous travaillerons sur de l'ADN, avec des sites supposés indépendants et un alignement dont les colonnes ont été vérifiées.

Introduction

La vraisemblance est la réponse probabiliste à une question simple : quelle est la probabilité d'observer cet alignement, étant donné un arbre, un modèle d'évolution et un jeu de paramètres ? Cette formulation renverse la perspective habituelle. Une distance résume une paire de séquences par un nombre ; une vraisemblance utilise toutes les colonnes de l'alignement, y compris celles qui ne varient pas, et fait dépendre chaque site de l'arbre entier. Elle permet de comparer des topologies, d'estimer des longueurs de branches et de choisir un modèle dans un cadre statistique unique.

L'enjeu biologique est de ne pas confondre une histoire plausible avec une histoire bien soutenue. Des lignées à évolution rapide, des compositions biaisées et des sites soumis à sélection produisent des convergences qui trompent les méthodes de distances. Un modèle explicite peut absorber une partie de ces effets, mais il introduit ses propres hypothèses : stationnarité, homogénéité, réversibilité, indépendance des sites. La vraisemblance ne supprime pas ces hypothèses ; elle les rend visibles et permet de les tester, au moins partiellement.

L'enjeu computationnel est triple. Il faut d'abord calculer une probabilité pour un site, puis multiplier sur des milliers de sites sans underflow. Il faut ensuite optimiser des longueurs de branches et des paramètres de modèle, souvent sous contraintes. Il faut enfin explorer un espace de topologies dont la taille explose : pour $n$ taxons, le nombre d'arbres binaires non enracinés vaut $(2n-5)!!$ , soit plus de deux millions dès dix taxons. Les algorithmes réels combinent donc élagage, optimisation numérique et heuristiques de recherche locale.

La sélection de modèle ajoute une couche. Un modèle plus riche augmente presque toujours la log-vraisemblance ; les critères d'information, comme le BIC de ModelFinder, pénalisent les paramètres. Les supports mesurent la répétabilité d'une bipartition sous rééchantillonnage, pas sa vérité. Cette section construit ces notions, puis les applique à un alignement synthétique de quatre taxons.

Théorie

Données, arbre, modèle et paramètres

Notons $D$ l'alignement, $T$ la topologie de l'arbre avec ses longueurs de branches, et $M$ le modèle d'évolution. Les paramètres $\theta$ regroupent les longueurs de branches, les échangeabilités et les fréquences du modèle, ainsi que les paramètres d'hétérogénéité des taux. La vraisemblance s'écrit :

\mathcal{L}(\theta) = P(D \mid T, M, \theta).

Cette quantité n'est pas la probabilité que l'arbre soit vrai. Elle ne contient ni prior, ni probabilité a posteriori : elle mesure la compatibilité des données avec une hypothèse. Un arbre faux peut avoir une vraisemblance élevée si le modèle est mal choisi. La formulation conditionnelle impose de comparer des valeurs calculées sur le même alignement et avec le même modèle.

Sous un modèle réversible, la racine n'affecte pas la vraisemblance : on travaille sur un arbre non enraciné. Les longueurs de branches s'expriment en substitutions par site, pas en temps, et un arbre de gènes n'est pas automatiquement un arbre d'espèces.

La vraisemblance d'un site

On modélise l'évolution le long d'une branche par une chaîne de Markov à temps continu. Pour un modèle à quatre états, la matrice de taux $Q$ vérifie, pour $i \neq j$ , $Q_{ij} = \pi_j r_{ij}$ , où $\pi_j$ est la fréquence stationnaire de l'état $j$ et $r_{ij}$ l'échangeabilité. La diagonale est choisie pour que chaque ligne somme à zéro, et l'échelle globale est fixée pour que le taux moyen vaille un. La probabilité de passer de l'état $i$ à l'état $j$ en une branche de longueur $t$ est l'élément $(i,j)$ de la matrice exponentielle :

P(t) = e^{Qt}.

Pour un site donné, on observe un état à chaque feuille. L'élagage de Felsenstein calcule la probabilité de ces observations en remontant l'arbre. Soit $L_u(x)$ la vraisemblance conditionnelle du sous-arbre sous le nœud $u$ , sachant que $u$ porte l'état $x$ . Pour une feuille observée dans l'état $y$ , $L_u(x)$ vaut 1 si $x=y$ et 0 sinon. Pour un nœud interne de children $c_1, \dots, c_k$ , on a :

L_u(x) = \prod_{c} \sum_{y} P_{xy}(t_c) L_c(y).

À la racine, on moyenne sur les états avec les fréquences stationnaires : $P(D_i) = \sum_x \pi_x L_{\mathrm{root}}(x)$ . Cette récursion évite d'énumérer toutes les histories possibles sur l'arbre. Sa complexité est de l'ordre de $O(n s^2)$ par site, où $n$ est le nombre de nœuds et $s$ le nombre d'états (4 pour l'ADN, 20 pour les protéines). Pour $L$ sites, une évaluation complète coûte $O(L n s^2)$ .

Les sites sont indépendants parce que le modèle le postule, non parce que la biologie le garantit. L'hypothèse est fragile pour des tiges appariées ou des régions sélectionnées. Les modèles à hétérogénéité ajoutent une mixture : proportion de sites invariables ou catégories gamma. La vraisemblance d'un site devient une somme de vraisemblances conditionnelles, plus coûteuse mais plus réaliste.

Log-vraisemblance et optimisation

Le produit de milliers de probabilités devient très petit et dépasse rapidement la précision des flottants. On travaille donc avec la log-vraisemblance :

\ell(\theta) = \ln \mathcal{L}(\theta) = \sum_{i=1}^{L} \ln P(D_i \mid T, M, \theta).

La log-vraisemblance est négative ou nulle, et une valeur plus grande (moins négative) indique un meilleur ajustement. Les logiciels utilisent des facteurs d'échelle par site pour éviter l'underflow, puis additionnent les logarithmes. Le maximum de vraisemblance est la valeur de $\theta$ qui maximise $\ell$ , pas la valeur qui maximise chaque site séparément.

L'optimisation porte sur des dizaines de variables contraintes : branches positives, fréquences de somme un, forme gamma positive. Les algorithmes combinent optimisations unidimensionnelles et méthodes quasi-Newton, en partant d'arbres variés. Une log-vraisemblance élevée n'est donc pas une preuve d'optimalité globale : c'est le meilleur point trouvé par une procédure donnée.

Exemple travaillé à la main

Prenons deux séquences et deux sites : le premier site est identique, le second diffère. Sous Jukes-Cantor, pour une distance totale $d$ entre les deux séquences, la probabilité qu'un site soit identique vaut :

P(\text{identique}) = \frac{1}{4} + \frac{3}{4}e^{-4d/3},

et la probabilité qu'il diffère vaut $1 - P(\text{identique})$ . Pour $d = 0{,}1$ , on obtient $P(\text{identique}) \approx 0{,}9064$ et $P(\text{different}) \approx 0{,}0936$ . La vraisemblance des deux sites vaut donc $0{,}9064 \times 0{,}0936 \approx 0{,}0849$ , soit $\ell \approx -2{,}466$ .

Pour $d = 0{,}5$ , les probabilités valent environ $0{,}669$ et $0{,}331$ , et $\ell \approx -1{,}509$ . La valeur a augmenté : le modèle préfère une distance plus grande, car un site différent est mieux expliqué par une divergence modérée. Pour $d = 0{,}824$ , on obtient $P(\text{identique}) \approx 0{,}5$ et $P(\text{different}) \approx 0{,}5$ , donc $\mathcal{L} = 0{,}25$ et $\ell \approx -1{,}386$ . C'est le maximum pour ces deux sites : la proportion observée de différences vaut $0{,}5$ , et la correction de Jukes-Cantor donne précisément $d \approx 0{,}824$ . La vraisemblance retrouve la distance corrigée, mais elle fournit en plus une valeur comparable entre arbres et entre modèles.

Recherche d'arbres : NNI et SPR

Pour un nombre de taxons donné, l'espace des topologies est immense. Une recherche exhaustive est hors de portée dès une dizaine de taxons. Les logiciels utilisent des heuristiques de recherche locale. La plus simple est le plus proche voisin (nearest neighbor interchange, NNI) : on choisit une branche interne, on échange l'un des deux sous-arbres d'un côté avec l'un des deux sous-arbres de l'autre. Pour $n$ taxons, il y a $n-3$ branches internes et deux échanges possibles par branche, soit $2(n-3)$ voisins.

Le subtree pruning and regrafting (SPR) est plus général : on coupe un sous-arbre, puis on le réattache sur une autre branche. Le voisinage est bien plus vaste, de l'ordre de $O(n^2)$ , et chaque voisin demande une nouvelle optimisation. Le tree bisection and reconnection (TBR) va encore plus loin. Un algorithme glouton évalue les voisins, garde le meilleur, puis recommence jusqu'à ce qu'aucune amélioration ne soit trouvée. Il s'arrête souvent dans un optimum local.

IQ-TREE et RAxML-NG multiplient les départs et les perturbations : IQ-TREE alterne NNI et perturbations aléatoires, RAxML-NG part d'arbres de parcimonie et applique des SPR. La topologie finale dépend de la stratégie, du départ et du budget de calcul. Rapporter l'arbre sans la commande rend l'analyse difficile à reproduire.

ModelFinder, AIC et BIC

ModelFinder, intégré à IQ-TREE, parcourt des modèles de substitution (JC, K2P, HKY, GTR, etc.) combinés à des fréquences empiriques (+F), une proportion de sites invariables (+I), des catégories gamma (+G) ou des taux libres (+R). Pour chaque candidat, il estime la vraisemblance et calcule un critère d'information. Le BIC s'écrit :

\mathrm{BIC} = k \ln n - 2\ell,

où $k$ est le nombre de paramètres libres et $n$ le nombre de sites de l'alignement, et non le nombre de taxons. L'AIC s'écrit $\mathrm{AIC} = 2k - 2\ell$ . Dans les deux cas, une valeur plus petite est préférée. Le BIC pénalise plus fortement les modèles complexes lorsque $n$ grandit ; l'AIC tend à retenir des modèles plus riches.

Reprenons les valeurs illustratives de notre alignement : JC a $\ell = -28{,}5$ avec $k=5$ , K2P a $\ell = -27{,}8$ avec $k=6$ . Pour $n=16$ , le BIC de JC vaut $5 \ln 16 + 57 \approx 70{,}86$ ; celui de K2P vaut $6 \ln 16 + 55{,}6 \approx 72{,}24$ . K2P améliore la log-vraisemblance de 0,7, mais son paramètre supplémentaire coûte $1{,}4$ en BIC et $1{,}4$ en AIC : le modèle simple reste préféré. Ce n'est pas une vérité biologique, seulement un arbitrage entre ajustement et complexité sur 16 sites.

ModelFinder rapporte souvent plusieurs modèles proches. Un écart de BIC inférieur à 2 est faible, entre 2 et 6 positif, au-delà de 10 fort : ce sont des repères, pas des lois. Dans la zone d'incertitude, vérifiez si les conclusions changent ; un modèle complexe peut aussi être instable sur un alignement court.

Test du rapport de vraisemblance

Pour deux modèles imbriqués, on peut comparer leurs maxima par le rapport de vraisemblance. Si le modèle 2 contient le modèle 1 et ajoute $d$ paramètres, la statistique vaut :

\mathrm{LR} = 2(\ell_2 - \ell_1),

et suit approximativement une loi du khi-deux à $d$ degrés de liberté sous l'hypothèse nulle. Une p-value supérieure à 0,05 signifie que l'amélioration n'est pas significative au seuil choisi ; elle ne prouve pas que le modèle simple soit vrai. Le test exige des modèles imbriqués, un même alignement et des conditions de régularité. Les paramètres à la frontière, comme une proportion de sites invariables proche de zéro, rendent l'approximation fragile ; le BIC est alors plus prudent.

Dans notre exemple, $\mathrm{LR} = 2(-27{,}8 + 28{,}5) = 1{,}4$ avec $d=1$ . La p-value vaut environ 0,24 : on ne rejette pas JC au profit de K2P. Le test demande si l'amélioration est compatible avec le hasard de l'échantillonnage ; le BIC, si elle justifie le coût d'un paramètre.

UFBoot et bootstrap

Le bootstrap non paramétrique de Felsenstein rééchantillonne les sites avec remise : on tire $L$ colonnes, on reconstruit un arbre, puis on compte la fréquence d'une bipartition. Un support de 95 % signifie qu'elle est retrouvée dans environ 95 % des alignements rééchantillonnés, pas qu'elle a 95 % de chances d'être vraie.

Le bootstrap standard est coûteux : mille réplicats valent mille recherches complètes. L'ultrafast bootstrap (UFBoot), amélioré par UFBoot2, approxime ces réplicats à partir des log-vraisemblances de sites déjà calculées. Il est proche du bootstrap standard pour un coût bien plus faible, mais peut être anti-conservateur sur des branches courtes ou sous un modèle mal spécifié.

Les seuils usuels sont 95 % pour un support fort, 70 à 95 % pour un support modéré, moins de 70 % pour un support faible. Ces seuils dépendent du nombre de taxons et de sites ; ils ne dispensent pas de lire l'arbre. Un support élevé sur une branche courte peut disparaître avec quelques sites supplémentaires. Un support faible en profondeur peut refléter une vraie incertitude ou une attraction des branches longues.

Implémentation

Les commandes ci-dessous supposent les outils installés et exécutées depuis la racine du dépôt. IQ-TREE écrit un rapport .iqtree, un arbre .treefile et un journal .log. L'option -m MFP active ModelFinder puis la recherche d'arbre ; -B 1000 demande 1000 réplicats UFBoot ; -T 4 limite le calcul à quatre threads. La commande exacte utilisée pour l'étude de cas est :

iqtree2 -s datasets/phylogeny/demo.fasta -m MFP -B 1000 -T 4

RAxML-NG utilise une autre syntaxe. La commande suivante enchaîne la recherche et le bootstrap avec un modèle fixé ; elle est illustrative et doit être adaptée au modèle choisi :

raxml-ng --all --msa datasets/phylogeny/demo.fasta --model GTR+G \
  --bs-trees 1000 --threads 4 --prefix demo_ml

Le code Python ci-dessous rassemble les quatre fonctions de la section. Il ne lit pas de fichier : le rapport est une chaîne, ce qui rend les tests reproductibles. Les fonctions valident les entrées et refusent de fabriquer une valeur manquante.

Parser un rapport, choisir un modèle et comparer

import math
import re

NUMBER = r"[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?"

def parse_iqtree_report(text):
  if not isinstance(text, str):
      raise TypeError("the report must be a string")
  if not text.strip():
      raise ValueError("the report must be a non-empty string")
  def search(pattern):
      match = re.search(pattern, text, flags=re.IGNORECASE | re.MULTILINE)
      return match.group(1) if match else None
  model = search(r"Best(?:-fit)? model(?:\s+according to BIC)?\s*:\s*(\S+)")
  if model is None:
      raise ValueError("the report does not contain a best-fit model")
  log_likelihood = search(r"Log-likelihood of the tree\s*:\s*(" + NUMBER + ")")
  if log_likelihood is None:
      raise ValueError("the report does not contain a tree log-likelihood")
  def optional_int(pattern):
      value = search(pattern)
      return int(value) if value is not None else None
  def optional_float(pattern):
      value = search(pattern)
      return float(value) if value is not None else None
  return {
      "best_model": model,
      "log_likelihood": float(log_likelihood),
      "n_sequences": optional_int(r"Input data\s*:\s*(\d+)\s+sequences"),
      "n_sites": optional_int(r"sequences?\s+with\s+(\d+)\s+columns"),
      "n_params": optional_int(r"Number of free parameters[^:]*:\s*(\d+)"),
      "aic": optional_float(r"\bAIC\s+score\s*:\s*(" + NUMBER + ")"),
      "bic": optional_float(r"\bBIC\s+score\s*:\s*(" + NUMBER + ")"),
      "tree_length": optional_float(r"Total tree length[^:]*:\s*(" + NUMBER + ")"),
  }

def best_model(candidates):
  if not candidates:
      raise ValueError("at least one candidate model is required")
  seen = set()
  best_name = None
  best_bic = None
  for candidate in candidates:
      name = candidate.get("name")
      if not isinstance(name, str) or not name:
          raise ValueError("each candidate must have a non-empty name")
      if name in seen:
          raise ValueError("candidate names must be unique")
      seen.add(name)
      bic = float(candidate["bic"])
      if not math.isfinite(bic):
          raise ValueError("BIC must be finite")
      if best_bic is None or bic < best_bic:
          best_name, best_bic = name, bic
  return best_name

def compare_models(models, n_sites):
  if not models:
      raise ValueError("at least one model is required")
  rows = []
  for model in models:
      name = model["name"]
      log_likelihood = float(model["log_likelihood"])
      k = int(model["n_params"])
      rows.append({
          "name": name,
          "log_likelihood": log_likelihood,
          "n_params": k,
          "aic": 2.0 * k - 2.0 * log_likelihood,
          "bic": k * math.log(n_sites) - 2.0 * log_likelihood,
      })
  best = min(row["bic"] for row in rows)
  for row in rows:
      row["delta_bic"] = row["bic"] - best
  return sorted(rows, key=lambda row: row["bic"])

def likelihood_ratio(log_likelihood_1, log_likelihood_2, df):
  if df < 1:
      raise ValueError("df must be a positive integer")
  statistic = 2.0 * (log_likelihood_2 - log_likelihood_1)
  if statistic < -1e-9:
      raise ValueError("the more complex model must have the higher log-likelihood")
  return {"lr": max(0.0, statistic), "df": df}

report = """IQ-TREE 2.3.6
Input data: 4 sequences with 16 columns, 12 distinct patterns
Best-fit model according to BIC: JC
Log-likelihood of the tree: -28.5000 (s.e. 3.2000)
Number of free parameters (#branches + #model parameters): 5
AIC score: 67.0000
BIC score: 70.8631
Total tree length (sum of branch lengths): 0.5500
"""

parsed = parse_iqtree_report(report)
print(parsed["best_model"], round(parsed["log_likelihood"], 4), round(parsed["bic"], 4))
candidates = [
  {"name": "JC", "bic": parsed["bic"]},
  {"name": "K2P", "bic": 72.2355},
  {"name": "GTR+F+G4", "bic": 84.2985},
]
print(best_model(candidates))
for row in compare_models([
  {"name": "JC", "log_likelihood": -28.5, "n_params": 5},
  {"name": "K2P", "log_likelihood": -27.8, "n_params": 6},
  {"name": "GTR+F+G4", "log_likelihood": -26.9, "n_params": 11},
], 16):
  print(row["name"], round(row["aic"], 4), round(row["bic"], 4), round(row["delta_bic"], 4))
print(likelihood_ratio(-28.5, -27.8, 1))

SortieCtrl+Entrée pour exécuter · 10 s max

Aucune sortie. Exécutez le code pour voir le résultat.

Le code s'exécute dans votre navigateur (Web Worker isolé). Il n'a pas accès au serveur ; l'accès réseau est désactivé pendant l'exécution. Le premier clic sur « Exécuter » télécharge Pyodide depuis jsDelivr (plusieurs mégaoctets, version 0.26.4 épinglée) ; rien n'est chargé avant ce clic.

La sortie vérifiée le 2026-11-01 pour ce rapport illustratif est :

JC -28.5 70.8631
JC
JC 67.0 70.8629 0.0
K2P 67.6 72.2355 1.3726
GTR+F+G4 75.8 84.2985 13.4355
{'lr': 1.4000000000000021, 'df': 1}

Trois points critiques. Le BIC utilise le nombre de sites, pas de taxons. compare_models ne calcule ni p-value ni poids : le test du rapport de vraisemblance reste séparé. Le parseur est tolérant sur les options mais strict sur le modèle et la log-vraisemblance.

Les versions changent le format : IQ-TREE 2.3.6 écrit Best-fit model according to BIC et Log-likelihood of the tree, d'autres outils non. Testez le parseur sur les rapports réellement produits et conservez le rapport brut.

Étude de cas

Nous utilisons l'alignement synthétique datasets/phylogeny/demo.fasta, créé pour ce cours. Il contient quatre séquences de 16 pb, taxon_A à taxon_D, sans gap ni ambiguïté. Il ne représente aucune espèce réelle et ne doit pas servir à une conclusion biologique. Sa taille permet de vérifier les différences à la main et d'observer la fragilité d'un arbre sur peu de sites.

La question est la suivante : quel modèle le BIC retient-il, quelle topologie est trouvée, et que valent les supports ? Nous lançons la commande exacte, puis nous lisons le rapport. La sortie ci-dessous est illustrative : IQ-TREE n'a pas été exécuté dans cet environnement, et les nombres sont choisis pour respecter les relations entre log-vraisemblance, AIC et BIC. Ils servent à lire le format, pas à établir un résultat biologique.

IQ-TREE 2.3.6
Input data: 4 sequences with 16 columns, 12 distinct patterns
Best-fit model according to BIC: JC
Log-likelihood of the tree: -28.5000 (s.e. 3.2000)
Number of free parameters (#branches + #model parameters): 5
Akaike information criterion (AIC) score: 67.0000
Bayesian information criterion (BIC) score: 70.8631
Total tree length (sum of branch lengths): 0.5500
Tree in Newick format: ((taxon_A:0.08,taxon_B:0.08)78:0.11,taxon_C:0.14,taxon_D:0.14);

Le modèle JC est retenu. K2P améliorerait la log-vraisemblance de 0,7, mais coûte un paramètre de plus ; sur 16 sites, le BIC et l'AIC le pénalisent. Le test du rapport de vraisemblance donne $\mathrm{LR} = 1{,}4$ avec un degré de liberté, soit une p-value d'environ 0,24 : l'amélioration n'est pas significative. Cette conclusion vaut pour ce jeu, pas en général.

L'arbre place A et B ensemble, avec un support UFBoot de 78. A et B ne diffèrent qu'en une position, A et C également ; B et C diffèrent en deux positions et D est plus éloigné. Mais l'alternative A-C n'est séparée que par un site : avec 16 colonnes, un seul site peut faire basculer la topologie. Le support de 78 est modéré, et 22 % des rééchantillonnages ne retrouvent pas la bipartition. Ce n'est ni une probabilité de vérité ni un argument de parenté réelle.

La longueur totale vaut 0,55 substitution par site : 0,08 pour A et B, 0,11 pour la branche interne, 0,14 pour C et D. Ce sont des paramètres sous JC, pas des temps. Un modèle avec hétérogénéité produirait d'autres longueurs, et un alignement plus long serait nécessaire pour distinguer les topologies. La prudence s'impose sur la taille, le modèle et le support.

Erreurs fréquentes

Quiz

Exercices

Exercice guidé5.4-e1

Énoncé

Implémenter parse_iqtree_report(text) et best_model(candidates). Le parseur extrait le modèle, la log-vraisemblance, les nombres de séquences, de sites, de patterns et de paramètres, l'AIC, le BIC, la longueur de l'arbre, la forme gamma et la proportion de sites invariables. Les champs absents valent None ; une chaîne vide ou un rapport sans modèle ou sans log-vraisemblance lève ValueError. best_model reçoit des dictionnaires {"name", "bic"}, vérifie l'unicité des noms et renvoie le BIC minimal, avec l'ordre d'entrée comme départage.

Comment vérifier votre travail

Les tests se trouvent dans exercises/module-5/5.4-maximum-vraisemblance/test_5_4.py. Depuis la racine du dépôt :

python -m pytest exercises/module-5/5.4-maximum-vraisemblance -q

Corrigé verrouillé

Le corrigé est chargé uniquement après confirmation. Il n'est pas présent dans la page avant cette action.

Exercice approfondi5.4-e2

Énoncé

Implémenter likelihood_ratio(log_likelihood_1, log_likelihood_2, df) et compare_models(models, n_sites=None). La première calcule 2 * (logL2 - logL1), refuse un logL2 inférieur de plus de 1e-9 et renvoie la statistique, la p-value par khi-deux, le seuil et la décision. La seconde calcule AIC = 2k - 2 logL et BIC = k ln(n) - 2 logL, trie par BIC croissant et ajoute delta_bic et un poids de somme 1. Le nombre de sites est fourni une fois ou dans chaque modèle, et doit être identique.

Comment vérifier votre travail

Les tests se trouvent dans exercises/module-5/5.4-maximum-vraisemblance/test_5_4.py. Depuis la racine du dépôt :

python -m pytest exercises/module-5/5.4-maximum-vraisemblance -q

Corrigé verrouillé

Le corrigé est chargé uniquement après confirmation. Il n'est pas présent dans la page avant cette action.

Mini-projet5.4-projet

Énoncé

Mini-projet : rapport consolidé. Implémenter iqtree_report(text) qui valide un rapport IQ-TREE, extrait le modèle, la log-vraisemblance, les nombres de séquences, de sites et de paramètres, l'AIC et le BIC, puis lit les supports dans Tree in Newick format: .... Si le BIC manque mais que n_sites et n_params existent, le recalculer. Les supports entre 0 et 1 deviennent des pourcentages ; ceux supérieurs à 100 sont ignorés. Ajouter des avertissements pour un BIC recalculé ou incohérent, un alignement de moins de 100 sites, moins de quatre séquences, l'absence de supports et les supports inférieurs à 70. Renvoyer aussi summary. Un rapport sans modèle, sans log-vraisemblance ou sans moyen de calculer le BIC lève ValueError.

Comment vérifier votre travail

Les tests se trouvent dans exercises/module-5/5.4-maximum-vraisemblance/projet/test_projet_5_4.py. Depuis la racine du dépôt :

python -m pytest exercises/module-5/5.4-maximum-vraisemblance/projet -q

Corrigé verrouillé

Le corrigé est chargé uniquement après confirmation. Il n'est pas présent dans la page avant cette action.

Résumé

La vraisemblance est $P(D \mid T, M, \theta)$ ; elle se calcule par élagage, site par site, puis se multiplie sous forme de log-vraisemblance.
Les paramètres comprennent les longueurs de branches, les échangeabilités, les fréquences et l'hétérogénéité des taux ; la racine est arbitraire sous un modèle réversible.
La recherche d'arbres est heuristique : NNI échange des sous-arbres voisins, SPR élague et réattache, et plusieurs départs limitent les optima locaux.
ModelFinder choisit un modèle par BIC ou AIC ; le BIC vaut $k \ln n - 2\ell$ avec $n$ le nombre de sites, et pénalise la complexité.
Le test du rapport de vraisemblance ne s'applique qu'à des modèles imbriqués ; il complète le BIC mais ne le remplace pas.
UFBoot estime la fréquence des bipartitions sous rééchantillonnage des sites ; un support de 95 % n'est pas une probabilité de vérité.
Une analyse de vraisemblance doit rapporter l'alignement, le modèle, la commande, les versions, les supports et les limites biologiques.

Références

Felsenstein, J. (1981). Evolutionary trees from DNA sequences: a maximum likelihood approach. Journal of Molecular Evolution, 17(6), 368–376.
Felsenstein, J. (1985). Confidence limits on phylogenies: an approach using the bootstrap. Evolution, 39(4), 783–791.
Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A. & Jermiin, L. S. (2017). ModelFinder: fast model selection for accurate phylogenetic estimates. Nature Methods, 14(6), 587–589.
Hoang, D. T., Chernomor, O., von Haeseler, A., Minh, B. Q. & Vinh, L. S. (2018). UFBoot2: improving the ultrafast bootstrap approximation. Molecular Biology and Evolution, 35(2), 518–522.
Kozlov, A. M., Darriba, D., Flouri, T., Morel, B. & Stamatakis, A. (2019). RAxML-NG: a fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference. Bioinformatics, 35(21), 4453–4455.
Minh, B. Q., Schmidt, H. A., Chernomor, O., Schrempf, D., Woodhams, M. D., von Haeseler, A. & Lanfear, R. (2020). IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Molecular Biology and Evolution, 37(5), 1530–1534.
Documentation IQ-TREE 2.3.6. http://www.iqtree.org/doc/ ; documentation RAxML-NG 1.2.0. https://github.com/amkozlov/raxml-ng/wiki