Fundamentos matemáticos de la coalescencia en Biología

Traballo Fin de Grao Fundamentos matemáticos de la coalescencia en Biología Mar Vázquez Rabuñal 2020/2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Fundamentos matemáticos de la coalescencia en Biología Mar Vázquez Rabuñal Julio, 2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Coñecemento: Estadística e investigación operativa Título: Fundamentos matemáticos de la coalescencia en Biología Breve descrición do contido La teoría de la coalescencia estudia cómo las variantes genéticas de una población pueden haberse originado a partir de un ancestro común. En el caso más simple, se supone que no hay recombinación, ni selección natural, ni estructura de la población, lo que significa que es igualmente probable que cada variante haya pasado de una generación a la siguiente. El objetivo de este trabajo es introducir al alumno en el estudio de las herramientas matemáticas de la teoría de la coalescencia, que son una colección de modelos estocásticos utilizados para generar predicciones sobre patrones de variación genética, y saber hacer inferencias a partir de muestras de datos genéticos. Recomendacións Outras observacións iii Índice general Resumen viii Introducción xi 1. Conceptos preliminares 1 1.1. Conceptosbiológicos............................... 1 1.2. Conceptos matemáticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1. Distribuciones de probabilidad útiles en este trabajo . . . . . . . . . 8 2. Coalescencia en el modelo de Wright-Fisher 15 2.1. Modelo de Wright-Fisher . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1.1. Número de descendientes de un gen en una generación . . . . . . . . 17 2.2. Fundamentos de la coalescencia . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.1. Coalescencia en tiempo discreto . . . . . . . . . . . . . . . . . . . . . 20 2.2.2. Coalescencia en tiempo continuo . . . . . . . . . . . . . . . . . . . . 23 2.2.3. Modelo alternativo: Modelo de Moran . . . . . . . . . . . . . . . . . 26 2.3. Medidas del tamaño de una genealogía . . . . . . . . . . . . . . . . . . . . . 27 3. Coalescencia en el modelo de Wright-Fisher con mutaciones 35 3.1. Modelo de Wright-Fisher con mutaciones . . . . . . . . . . . . . . . . . . . . 35 3.1.1. Modelo de sitios infinitos . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2. Mutaciones y coalescencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.3. Algoritmos de generación de genealogías . . . . . . . . . . . . . . . . . . . . 40 3.4. Medidas de polimorfismos en una secuencia de ADN . . . . . . . . . . . . . 43 4. Aplicaciones de la teoría de la coalescencia 49 A. Códigos de R 57 A.1.Algoritmo1.................................... 57 v A.2.Algoritmo2.................................... 58 A.3.Algoritmo3.................................... 60 A.4. Diagrama de barras 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 A.5.Neandertales ................................... 62 Bibliografía 67 vi Figura 1.1: Componentes principales de un nucleótido. La mayor parte del ADN presente en los seres vivos se encuentra compactificado en unas estructuras con forma de X denominadas cromosomas. Según el número de cromosomas que presente un individuo se dirá que es haploide odiploide. Los individuos haploides son aquellos que presentan un único juego de cromosomas y que para reproducirse se duplican y dividen. Son haploides, por ejemplo, algunas bacterias. Los individuos diploides son aquellos que presentan dos juegos de cromosomas, uno procedente del padre y otro de la madre. El ser humano es diploide y presenta 23 pares de cromosomas, un par de ellos son los denominados cromosomas sexuales X e Y (X hace referencia al sexo femenino e Y al masculino). Si un ser humano presenta estos dos cromosomas como XX será una mujer y si presenta XY será un hombre. Ahora que ya se conoce cómo es la estructura del material genético a nivel orgánico, es importante centrarse en los conceptos con los que trata la genética de poblaciones. El primero de ellos, y el más importante, es el concepto de gen. Una definición estricta establece que un gen es la unidad funcional de la herencia que controla cada carácter de los seres vivos. A nivel estructural un gen no es más que una secuencia o segmento de ADN que codifica una determinada información. Al conjunto de todos los genes de una especie se le denomina genoma. Cada gen que lleva la información de un determinado carácter puede manifestarse de varias formas. A cada una de estas formas se les conoce como alelos. Un ejemplo clásico para explicar los alelos es el del color de los guisantes. El gen que lleva la información del color presenta dos alelos: uno que determina el color verde y otro el color amarillo. Según cuál de los alelos presente el gen, el guisante será verde o será amarillo. En ocasiones las secuencias de ADN experimentan mutaciones, que no son más que cambios en la secuencia de nucleótidos que lo conforman. Este tipo de modificaciones 2 pueden tener efectos muy diversos: desde que sean imperceptibles, hasta que supongan una ventaja o una desventaja para el organismo que la presenta. Dentro de la genética poblacional es muy común el estudio de la historia de una determinada población, la llamada genealogía, que no es más que el análisis de las distintas generaciones de una población. Para facilitar este estudio se suele emplear la estructura de árbol genealógico que conecta los parentescos de los individuos, genes o secuencias de ADN que se estén considerando. Una palabra muy común dentro de este contexto es la de linaje, que hace referencia a la línea de antepasados de un determinado individuo. En la teoría de la coalescencia nos interesa conocer el ancestro común más reciente de una determinada muestra de genes. Dibujando el árbol genealógico sobre los genes de estudio, podremos visualizar con facilidad dónde se encuentra este ancestro común, aunque después nos será necesario establecer un modelo matemático para poder analizar todo el proceso más a fondo. Teniendo en cuenta todos estos conceptos propios de la estructura biológica del material genético, así como las ideas básicas de la genética poblacional, contamos con los conocimientos necesarios, a nivel biológico, para adentrarnos en la teoría de la coalescencia. Pese a todo, vamos a hacer una pequeña presentación de conceptos matemáticos que nos serán de utilidad para modelizar esta teoría y entenderla. Para ello nos basaremos principalmente en el libro de Sheldon Ross: A First Course in Probability ([1]). 1.2. Conceptos matemáticos Consideremos un experimento aleatorio, es decir, un experimento del que no se puede saber con certeza cuál va a ser su resultado. Al conjunto de todos los posibles resultados de un experimento se le conoce como espacio muestral del experimento y se denota por Ω. A cualquier subconjunto del espacio muestral se le denomina evento osuceso (E), es decir, un evento es un conjunto de posibles resultados de un experimento. Podemos definir la probabilidad de un evento Edel espacio muestral Ω,P(E), de forma que satisface los tres axiomas de Kolmogorov. Axioma 1: 0≤P(E)≤1 Axioma 2: P(Ω) = 1 Axioma 3: Para cualquier secuencia de eventos mutuamente excluyentes E1, E2,· · · (es decir, eventos para los cuales Ei∩Ej=∅cuando i6=j), se tiene: P(E1∪E2∪ · · · ) = X i P(Ei) En muchas ocasiones, a la hora de llevar a cabo un experimento, no nos interesan los 3 propios resultados del mismo, sino cantidades definidas a partir de ellos. Estas cantidades, que formalmente son funciones reales sobre el espacio muestral, se denominan variables aleatorias. Una variable aleatoria es una función Xde la forma: X: Ω −→ R. Como el valor de una variable aleatoria viene determinado por el resultado del experimento, podemos asignar probabilidades a los posibles valores de la variable aleatoria. Para una variable aleatoria Xse define su función de distribución Fde la siguiente manera: F(x) = P(X≤x),con − ∞ <x<∞(1.1) Por lo tanto, la función de distribución es simplemente aquella función que determina para todos los valores reales de x, la probabilidad de que la variable aleatoria sea menor o igual que x. Una variable aleatoria que presenta un número finito o infinito numerable de valores posibles se dice que es discreta. Dada una variable discreta X, se define la función de masa de probabilidad p(x)de Xcomo: p(x) = P(X=x)(1.2) Por otro lado, diremos que una variable aleatoria Xes continua cuando puede tomar cualquier valor en algún intervalo (o intervalos) del conjunto de los número reales. En esta situación se denomina función de densidad de probabilidad de Xa una función no negativa f, definida para todos los números reales x∈(−∞,∞), que satisface que para cualquier conjunto Bde números reales se tiene: P(X∈B) = ZB f(x)dx (1.3) Además, esta función debe cumplir: P[X∈(−∞,∞)] = Z∞ −∞ f(x)dx = 1 (1.4) Un concepto muy importante en la teoría de la probabilidad es el de esperanza de una variable aleatoria. Si Xes una variable aleatoria discreta con función de masa de probabilidad p(x), entonces su valor esperado o esperanza E(X)viene dado por: E[X] = X x xp(x)(1.5) Por otro lado, si Xes una variable aleatoria continua se define la esperanza de Xcomo: E[X] = Z∞ −∞ xf(x)dx (1.6) 4 Una propiedad muy empleada de la esperanza es que, dada una función real g, si tenemos una variable discreta Xo una variable continua Y, se verifican: E[g(X)] = X x g(x)p(x)E[g(Y)] = Z∞ −∞ g(y)f(y)dy (1.7) Podemos ver la demostración de esta propiedad en el libro A First Course in Probability de Sheldon Ross ([1]). Usando esta última propiedad vemos que, dadas aybdos constantes determinadas, se tiene que: E[a+bX] = a+bE[X]. La esperanza de una variable aleatoria no nos informa sobre la dispersión de los posibles valores de X. Una medida que sí que nos aporta esta información es la llamada varianza. La varianza de una variable aleatoria Xse define de la siguiente manera: V ar(X) = E(X−µ)2(1.8) donde µ=E[X]. Una forma alternativa de escribir V ar(X)es: V ar(X) = EX2−(E[X])2(1.9) Veamos que esta forma es equivalente a la dada en la definición, utilizando las propiedades de la esperanza: V ar(X) = E(X−µ)2=EX2−2µX +µ2 =EX2−2µE[X] + µ2=EX2−µ2(1.10) Otra medida interesante es la llamada covarianza de dos variables aleatorias XeY, que nos da información sobre la relación existente entre ellas. Se define como: Cov(X, Y ) = E[(X−E[X])(Y−E[Y])] = E[XY ]−E[X]E[Y](1.11) Para llegar a la segunda igualdad se emplean las propiedades básicas de la esperanza de una variable aleatoria. Algunas propiedades inmediatas de esta cantidad son que Cov(X, Y ) = Cov(Y, X)y que V ar(X) = Cov(X, X). Otra propiedad sencilla es la que sigue: Cov   n X i=1 Xi, m X j=1 Yj = n X i=1 m X j=1 Cov(Xi, Yj)(1.12) 5 Probémosla escribiendo para abreviar E[Xi] = µiyE[Yj] = νj. Además, tengamos en cuenta que E[Pn i=1 Xi] = Pn i=1 µiy que EhPm j=1 Yji=Pm j=1 νj. Así: Cov   n X i=1 Xi, m X j=1 Yj =E  n X i=1 Xi− n X i=1 µi!  m X j=1 Yj− m X j=1 νj   = n X i=1 m X j=1 E[(Xi−µi)(Yj−νj)] = n X i=1 m X j=1 Cov(Xi, Yj)(1.13) Ahora que sabemos esto podemos demostrar la siguiente propiedad: V ar n X i=1 Xi!= n X i=1 V ar(Xi)+2XX i<j Cov(Xi, Xj)(1.14) Tenemos: V ar n X i=1 Xi!=Cov   n X i=1 Xi, n X j=1 Xj = n X i=1 n X j=1 Cov (Xi, Xj) = n X i=1 V ar(Xi)+2XX i<j Cov(Xi, Xj)(1.15) En el último paso hemos usado que cada par de índices i, j, con i6=j, aparece dos veces en la doble suma anterior. Definamos ahora un concepto también esencial que es el de independencia de variables aleatorias. Dadas dos variables aleatorias XeY, decimos que son independientes si para cualquier par de conjuntos de números reales AyBse tiene que: P(X∈A, Y ∈B) = P(X∈A)P(Y∈B)(1.16) Consideremos la función de densidad conjunta de dos variables continuas XeY,f(x, y), definida por: P[(X, Y )∈C] = Z Z (x,y)∈C f(x, y)dxdy (1.17) donde C⊂R2. Se tiene que XeYson independientes si, y solo si, f(x, y) = fX(x)fY(y)para todo x, y ∈R, donde fXyfYson las funciones de densidad marginales de Xy de Y, respectivamente. 6 En el caso de que XeYsean variables discretas, la condición de independencia es equivalente a p(x, y) = pX(x)pY(y), para todo x, y ∈R. En esta situación, es decir, si XeYson independientes se tiene que Cov(X, Y )=0. Veámoslo probando que E[XY ] = E[X]E[Y]y considerando que XeYson continuas (si fuesen discretas el razonamiento sería muy parecido): E[XY ] = Z∞ −∞ Z∞ −∞ xyf(x, y)dxdy =Z∞ −∞ Z∞ −∞ xyfX(x)fY(y)dxdy =Z∞ −∞ xfX(x)dx Z∞ −∞ yfY(y)dy =E[X]E[Y](1.18) Por lo tanto, si X1, X2,· · · , Xnson variables aleatorias mutuamente independientes, entonces, aplicando que la covarianza entre dos de ellas es siempre 0, se tiene, observando la ecuación (1.15), que: V ar n X i=1 Xi!= n X i=1 V ar(Xi)(1.19) Un concepto que también usaremos a lo largo del trabajo es el de probabilidad condicionada. Sean dos eventos AyBdel espacio muestral Ωcon P(A)>0, entonces la probabilidad condicionada de Bdado A, es definida como: P(B|A) = P(A∩B) P(A)(1.20) En relación a esto tenemos el teorema de las probabilidades totales, que usaremos más adelante, y que nos dice: dados A1, A2, . . . , Akformando una partición de Ω(es decir, Ai∩Aj=∅si i6=jyPk i=1 Ai= Ω) y teniendo 0< P(Ai)<1para todo i, entonces dado un evento Ben Ω: P(B) = k X i=1 P(B|Ai)P(Ai)(1.21) Algo que también nos va a ser de utilidad es calcular la función de distribución de la suma de dos variables aleatorias. Consideremos así dos variables aleatorias XeYy calculemos la función de distribución de Z=X+Y, es decir, la convolución de las funciones de distribución de XeY: FZ(z) = P(X+Y≤z) = Z∞ −∞ P(X+Y≤z|Y=y)fY(y)dy =Z∞ −∞ P(X≤z−y|Y=y)fY(y)dy =Z∞ −∞ FX|Y(z−y)fY(y)dy (1.22) donde FX|Yes la función de distribución condicionada de Xdado Y=y. Si XeYson independientes, entonces FX|Y=FX. De esta forma: FZ(z) = Z∞ −∞ FX(z−y)fY(y)dy (1.23) 7 Derivando esta expresión respecto a zllegamos a la función de densidad de la suma: fZ(z) = Z∞ −∞ fX(z−y)fY(y)dy (1.24) Ahora que ya conocemos todos estos conceptos y razonamientos esenciales de la teoría de la probabilidad, vamos a presentar algunas de las distribuciones de probabilidad de variables aleatorias más comunes y que nos van a ser de gran utilidad más adelante. 1.2.1. Distribuciones de probabilidad útiles en este trabajo Distribución de Bernoulli Sea X una variable aleatoria que puede tomar los valores 0y1, haciendo referencia al fracaso o al éxito, respectivamente. La probabilidad de tener un éxito es py la de obtener un fracaso es 1−p. Decimos que una variable aleatoria Xsigue una distribución de Bernoulli, X∼Ber(p), si su función de masa de probabilidad viene dada por: P(X=x) = px(1 −p)1−x,con x= 0,1(1.25) Podemos calcular fácilmente su esperanza y su varianza: E[X] = X x∈{0,1} x·P(X=x) = 0 ·P(X= 0) + 1 ·P(X= 1) = p(1.26) V ar(X) = E[X2]−E[X]2=X x∈{0,1} x2·P(X=x)−p2=p−p2=p(1 −p)(1.27) Distribución binomial Consideremos ahora nintentos independientes de Bernoulli de parámetro p, es decir, cada uno de ellos tiene probabilidad de éxito py probabilidad de fracaso 1−p. Si X representa el número de éxitos que tienen lugar en los nintentos, entonces Xse dice que sigue una distribución de probabilidad binomial con parámetros (n, p),X∼Bi(n, p). La función de masa de probabilidad de una variable aleatoria binomial de parámetros (n, p)viene dada por: P(X=x) = n xpx(1 −p)n−x,con x= 0,1,2, . . . , n (1.28) 8 Obtengamos la esperanza y la varianza de una variable aleatoria de este tipo. Para ello calculemos antes de nada E[Xk]con k∈ {1,2,3, . . .}. E[Xk] = n X x=0 xkn xpx(1 −p)n−x= n X x=1 xkn xpx(1 −p)n−x =np n X x=1 xk−1n−1 x−1px−1(1 −p)n−x =np n−1 X i=0 (i+ 1)k−1n−1 ipi(1 −p)n−i−1 =npE[(Y+ 1)k−1](1.29) donde hemos usado que i=x−1y donde Yes una variable aleatoria binomial de parámetros (n−1, p). Sustituyendo para k= 1 llegamos a que E[X] = np y, por tanto, para k= 2 tenemos: E[X2] = npE[Y+ 1] = np [(n−1)p+ 1] (1.30) Así, la varianza viene dada por: V ar(X) = E[X2]−E[X]2=np [(n−1)p+ 1] −(np)2=np(1 −p)(1.31) Distribución geométrica Consideremos una serie de intentos independientes, cada uno de ellos con dos posibles resultados: éxito (con probabilidad p) o fracaso (con probabilidad q= 1 −p). Sea X la variable que describe el número de intentos necesarios para conseguir el primer éxito, decimos que Xsigue una distribución geométrica de parámetro p,X∼Geo(p). Esta variable tiene una función de masa de probabilidad dada por: P(X=x) = p(1 −p)x−1,con x= 1,2, . . . (1.32) Para las variables que siguen esta distribución podemos calcular la esperanza como se muestra a continuación: E[X] = ∞ X x=1 xpqx−1= ∞ X x=1 (x−1 + 1)pqx−1= ∞ X x=1 (x−1)pqx−1+ ∞ X x=1 pqx−1 =q ∞ X i=0 ipqi−1+ 1 = qE[X]+1 (1.33) 9 Despejando esto llegamos a que E[X] = 1 p. Calculemos ahora E[X2]para poder obtener después V ar(X). E[X2] = ∞ X x=1 x2pqx−1= ∞ X x=1 (x−1 + 1)2pqx−1 = ∞ X x=1 (x−1)2pqx−1+ ∞ X x=1 2(x−1)pqx−1+ ∞ X x=1 pqx−1 =q ∞ X i=0 i2pqi−1+ 2q ∞ X i=0 ipqi−1+ 1 =qE[X2]+2qE[X] + 1 = qE[X2] + 2q p+ 1 (1.34) Despejando de esta última expresión llegamos a que E[X2] = q+1 p2. Obtengamos ahora V ar(X): V ar(X) = E[X2]−E[X]2=q+ 1 p2−1 p2=1−p p2(1.35) Para probar una serie de propiedades interesantes, nos va a ser útil conocer P(X≤x0): P(X≤x0) = x0 X x=1 pqx−1=p x0−1 X i=0 qi=p1−qx0 1−q= 1 −(1 −p)x0(1.36) Ahora que ya sabemos esto, vamos a considerar la propiedad que hace referencia a la falta de memoria de esta distribución: P(X > x2+x1|X > x1) = P(X > x2)(1.37) Para entender bien esta propiedad de manera general imaginémonos que Xhace referencia al número de veces que sale cara antes de que salga cruz al tirar una moneda. La probabilidad de que salga cruz cuando han salido x1+x2caras, sabiendo que ya han salido x1caras, es la misma que la probabilidad inicial de que salga cruz después de x2caras. Es decir, se ha “olvidado” que ya han salido x1caras. Veamos su demostración: P(X > x2+x1|X > x1) = P(X > x2+x1, X > x1) P(X > x1)=P(X > x2+x1) P(X > x1) =(1 −p)x2+x1 (1 −p)x1= (1 −p)x2=P(X > x2)(1.38) Distribución de Poisson Sea Xuna variable aleatoria que mide el número de eventos que tienen lugar en un determinado intervalo. Se dice que Xsigue una distribución de Poisson de parámetro λ, 10 X∼Po(λ), si λes el número promedio de veces que se espera que ocurra el fenómeno en un intervalo dado. Se tiene que la función de masa de probabilidad es: P(X=x) = e−λλx x!,con x= 0,1,2, . . . (1.39) donde xes el número veces que ocurre el fenómeno. Podemos obtener la esperanza y la varianza de una variable aleatoria que siga esta distribución: E[X] = ∞ X x=0 xe−λλx x!=e−λ ∞ X x=1 λx (x−1)! =λe−λ ∞ X j=0 λj j!=λe−λeλ=λ(1.40) Para determinar la varianza demos primero el valor de E[X2]: E[X2] = ∞ X x=0 x2e−λλx x!=λ ∞ X x=1 xe−λλx−1 (x−1)! =λ ∞ X i=0 (i+ 1)e−λλi i! =λ"∞ X i=0 ie−λλi i!+ ∞ X i=0 e−λλi i!#=λ(λ+ 1) (1.41) De esta forma, restando, tenemos: V ar(X) = λ Distribución exponencial Sea Xuna variable continua que sigue una distribución exponencial de parámetro λ, X∼Exp(λ). Su función de densidad viene dada por: f(x) = (λe−λx si x ≥0 0si x < 0(1.42) La función de distribución se obtiene de la siguiente manera: F(x) = P(X≤x) = Zx 0 λe−λxdx = 1 −e−λx,con x≥0(1.43) Podemos obtener la esperanza y la varianza de una variable aleatoria que sigue una distribución exponencial. Para ello calculemos antes de nada E[Xk], con k > 0. E[Xk] = Z∞ 0 xkλe−λxdx =−xke−λx ∞ 0 +Z∞ 0 kxk−1e−λxdx = 0 + k λE[Xk−1] = k λE[Xk−1](1.44) Sustituyendo para k= 1 y después para k= 2 llegamos a: E[X] = 1 λE[X2] = 2 λ2(1.45) 11 V ar(Vi)=2N·1 2N·1−1 2N= 1 −1 2N(2.3) Que la media de Visea 1es consecuencia de que el tamaño de la población es constante: si la media del número de descendientes de un gen fuese mayor o menor que 1, entonces indicaría que la población está aumentando o disminuyendo de tamaño, respectivamente. Además, la distribución conjunta del número de descendientes de los 2Ngenes en una generación dada, es una multinomial de parámetros 2Nyp1=p2=· · · =p2N=1 2N. Los posibles sucesos característicos de esta distribución son que un determinado gen i(con 1≤i≤2N) sea el padre de un miembro de la siguiente generación. Así, podemos calcular la probabilidad de que el gen 1tenga v1descendientes en la generación t+ 1, el 2tenga v2y así sucesivamente hasta que el gen 2Ntenga v2Ndescendientes en esa generación (recordemos que se tiene que verificar que v1+v2+· · · +v2N= 2N). Esta probabilidad viene dada por: P(V1=v1,· · · , V2N=v2N) = (2N)! v1!· · · v2N!1 2Nv1 · · · 1 2Nv2N =(2N)! v1!· · · v2N!1 2N2N (2.4) También podemos obtener la covarianza y el coeficiente de correlación del número de descendientes de dos genes iyj: Cov(Vi, Vj) = −2Npipj=−1 2N(2.5) Cor(Vi, Vj) = Cov(Vi, Vj) pV ar(Vi)V ar(Vj)=−1 2N−1(2.6) Vemos que para valores grandes de 2N,ViyVjpresentan una correlación baja. Además, es de esperar que la covarianza entre ViyVjsea negativa pues, si el gen ideja muchos descendientes en la nueva generación, el gen jes más probable que deje pocos. Esto es una consecuencia de que el tamaño de la población sea constante e igual a 2N. Cuando el número 2Nes muy elevado, podemos ver que la distribución de Vitiende a una Poisson con media y varianza 1, Po(1): l´ım 2N→∞ P(Vi=v) = l´ım 2N→∞ 2N v 1 2Nv1−1 2N2N−v =1 v!l´ım 2N→∞ 2N! (2N−v)! 1 (2N−1)v1−1 2N2N =1 v!e−1(2.7) 18 De esta forma, podemos expresar de manera informal que cuando 2Nes suficientemente grande, la distribución que sigue el número de descendientes de un gen ien la generación siguiente es aproximadamente una Poisson de parámetro 1: P(Vi=v)≈1 v!e−1(2.8) Teniendo esto en cuenta, vemos que la probabilidad de que un gen no deje descendientes en una determinada generación es de aproximadamente e−1≈0,37 y, por tanto, hay una probabilidad de 0,63 de que sí que los deje. Dada, por ejemplo, una población de tamaño 2N= 10000, si estudiamos el número de genes ancestrales hace t= 15 generaciones veremos que tan solo tenemos 0,631510000 ≈10 genes. Es decir, el resto de genes, unos 9990, han perdido su linaje en estas 15 generaciones. Estos razonamientos solamente son válidos si estamos considerando un número muy grande de genes, en otro caso, si buscamos el ancestro común de un grupo pequeño de todos estos genes recurriríamos a la teoría de la coalescencia. Para ilustrar esto consideremos una población de 6genes de los cuales buscamos el ancestro común más reciente (MRCA) del primero, del tercero y del sexto. Como se observa en la Figura 2.2, dos generaciones hacia atrás el primer gen y el tercero encuentran su ancestro común y cuatro generaciones hacia atrás lo encuentran los tres genes que estábamos considerando. Presente Pasado Figura 2.2: Genealogía de tres genes en una población que presenta en total seis. Modelizando matemáticamente las distintas variables que intervienen en estos procesos, podremos llegar a conocer cómo es la estructura general de las genealogías de este tipo y algunas cantidades que resultan de interés, como el número de linajes que hay en un determinado instante de tiempo. 19 2.2. Fundamentos de la coalescencia Como ya hemos mencionado, es interesante estudiar cuánto tiempo pasa hasta que un conjunto de genes encuentran su ancestro común más reciente. Si medimos el tiempo en forma discreta, es decir, en generaciones, y consideramos la variable aleatoria Tque da información sobre el tiempo de espera hasta que aparece el MRCA de dos genes tenemos que Tsigue una distribución geométrica. El éxito tendrá lugar si se encuentra el ancestro común de esos dos genes, y el fracaso si no es así (con probabilidades py1−prespectivamente). Así, T∼Geo(p): P(T=t) = p(1 −p)t−1con t= 1,2, . . . (2.9) Que T=timplica que ha habido t−1fracasos (t−1generaciones en las que no se ha encontrado un ancestro común) antes del primer éxito, el encuentro de su MRCA. 2.2.1. Coalescencia en tiempo discreto Una vez que ya conocemos qué distribución sigue la variable que hace referencia al número de generaciones que pasan hasta que se encuentra el ancestro común más reciente de dos genes, podemos estudiar cómo es este proceso de coalescencia en muestras de dos genes y después generalizarlo para muestras con un número mayor de genes. Coalescencia para una muestra de dos genes Lo primero que haremos será estudiar la distribución del tiempo de espera hasta que se obtenga el MRCA de dos genes en un modelo haploide con 2Ngenes. La probabilidad de que estos dos genes encuentren un ancestro común en la primera generación hacia atrás en el tiempo es de 1 2Npues el primer gen puede elegir a su padre con libertad, pero el segundo debe escoger el mismo que el primero, es decir, solo puede escoger 1de las 2Nposibilidades existentes. La probabilidad de que dos genes tengan distintos ancestros es, por tanto, 1−1 2N. Como las selecciones de los padres en distintas generaciones son procesos independientes, la probabilidad de que dos genes encuentren un ancestro común tgeneraciones atrás en el tiempo es: 1−1 2Nt−11 2N(2.10) Esto representa que en las primeras t−1generaciones no se encuentra el ancestro común pero en la t-ésima sí. Por lo tanto, el tiempo de coalescencia hasta que dos genes encuentren su MRCA, que denotaremos T2, sigue una distribución geométrica con parámetro 1 2N: P(T2=t) = 1−1 2Nt−11 2Ncon t= 1,2, . . . (2.11) 20 El valor esperado de T2viene dado por: E(T2) = 1 1 2N = 2N. Esto indica que si nos vamos tantas generaciones hacia atrás en el tiempo, como genes hay en la población, se espera encontrar el ancestro común más reciente de dos genes. Coalescencia en una muestra de n genes Generalicemos ahora lo obtenido en el apartado anterior para una muestra de ngenes, donde nse considera mucho menor que 2N, el tamaño de la población. Empecemos obteniendo la distribución del tiempo de espera hasta que k(≤n)genes presenten menos de k linajes ancestrales, es decir, hasta que alguno de los kgenes comparta un ancestro común con otro gen de ese conjunto en la anterior generación. La probabilidad de que kgenes tengan kancestros distintos en la anterior generación se obtiene de forma similar al caso de dos genes: el primer gen puede escoger libremente entre los 2Ngenes, el segundo tiene que escoger un padre distinto y entonces solo puede escoger entre 2N−1, el tercero entre 2N−2y así se sigue hasta que el último puede escoger entre 2N−(k−1). Por tanto, esta probabilidad viene dada por: (2N−1) 2N (2N−2) 2N· · · (2N−k+ 1) 2N=1−1 2N1−2 2N· · · 1−k−1 2N = 1 − k−1 X i=1 i 2N+O1 N2 = 1 −k 21 2N+O1 N2(2.12) Para pasar de la primera línea a la segunda hemos realizado los productos y agrupado los sumandos con potencias mayores o iguales a 1 N2en el término O1 N2, verificándose entonces que l´ımN→∞ O1 N2 1 N2 =cte. Para pasar de la segunda línea a la tercera hemos usado que Pk−1 i=1 i=k 2=k(k−1) 2. Probemos esta igualdad por inducción en k, siendo k≥2. Paso 1: Estudiemos la validez de la fórmula para el caso k= 2: 2(2 −1) 2= 1 (2.13) Paso 2: Consideremos que la propiedad se cumple para un determinado ky veamos que se verifica también para k+ 1. Así, la hipótesis de inducción será: k−1 X i=1 i=k(k−1) 2(2.14) 21 Veamos qué ocurre para k+ 1: k X i=1 i= k−1 X i=1 i+k=k(k−1) 2+k=k(k−1) + 2k 2=(k+ 1)((k+ 1) −1) 2(2.15) Con esto concluimos la prueba y establecemos que la fórmula es válida para todo número entero kmayor o igual que 2. Como se asume que nes mucho menor que N, los términos con potencias de 1/N2o mayores, es decir, O1 N2, son despreciables y pueden ser ignorados. Esta aproximación es equivalente a ignorar la posibilidad de que más de un par de genes encuentren un ancestro común en la misma generación. Por lo tanto, cuando nes mucho menor que N, la probabilidad de que no se produzca ninguna coalescencia, es decir, que dados kgenes todos tengan un ancestro distinto en la anterior generación, es: 1−k 21 2N(2.16) y, por tanto, la probabilidad de que tenga lugar una coalescencia en una generación dada es: k 21 2N(2.17) En consecuencia, denotando por Tkel tiempo de coalescencia hasta que dos genes, de los kconsiderados, encuentren su MRCA, se tiene que la probabilidad de que dos de los kgenes considerados encuentren un ancestro común tgeneraciones en el pasado (con t= 1,2, . . . ) es, aproximadamente: P(Tk=t)≈1−k 21 2Nt−1k 21 2N(2.18) Por lo tanto, vemos que la variable aleatoria Tksigue, aproximadamente, una distribución geométrica con parámetro (k 2) 2N. En la Figura 2.3 hemos considerado 5genes y vemos el árbol de coalescencia asociado a su genealogía. Con este árbol buscamos comprender bien a qué nos referimos cuando hablamos de Tk, es decir, del tiempo en el que hay kancestros de los 5genes seleccionados. En la imagen vemos diferenciadas las distintas etapas Tkcon k= 2,3,4,5. 22 1 2 4 3 5 T5 T4 T3 T2 Figura 2.3: Árbol de coalescencia con los distintos Tkseñalados, siendo k= 2,3,4,5. La exactitud de las aproximaciones realizadas para poblaciones de gran tamaño, lleva a una formulación de la coalescencia en la que se emplea un modelo basado en la continuidad del tiempo y que es independiente de 2N. Es aquí donde introduciremos que el tiempo (continuo) hasta encontrar el MRCA de dos genes sigue una distribución exponencial. 2.2.2. Coalescencia en tiempo continuo Una de las formas más naturales de introducir el tiempo en escala continua en la teoría de la coalescencia es considerando que una unidad de tiempo se corresponde con el tiempo medio que tardan dos genes en encontrar un ancestro común (ya hemos visto que son 2Ngeneraciones). Esta transformación del tiempo hace que la coalescencia se vuelva independiente del tamaño de la población. Para derivar el proceso continuo de la coalescencia, se considera tc=t 2Ndonde t es el tiempo medido en generaciones. Obviamente, se puede pasar el tiempo continuo a generaciones simplemente despejando t; si el valor que se obtiene t= 2Ntcno es un entero entonces tse trunca al entero menor más cercano. Ahora veremos que la distribución geométrica, empleada en el apartado anterior, puede 23 ser aproximada por una distribución exponencial. Sea T∼Geo(p), se tiene: P(T > t) = 1 −P(T≤t) = 1 −1−(1 −p)t= (1 −p)t(2.19) Considerando a=p2N, podemos reescribir (1 −p)tcomo: 1−p2N 2N2Nt 2N =1−a 2Ntc2N(2.20) Así, en el límite de 2Nmuy grande se tiene: l´ım 2N→∞ P(T > t) = l´ım 2N→∞ PT 2N> tc= l´ım 2N→∞ 1−a 2Ntc2N=e−atc(2.21) De esta forma, la variable que definiremos por Tc=T 2N, sigue una distribución exponencial con parámetro aen el límite de 2Ntendiendo a ∞. Denotaremos el tiempo hasta que kgenes tengan k−1ancestros en el caso continuo como Tc k, que es simplemente una variable dada por Tk 2N, donde Tkes el tiempo medido en generaciones, para que kgenes tengan k−1ancestros. Recordemos que Tk∼Geo (k 2) 2N y, por lo que acabamos de ver, se tiene que: Tc k∼Exp k 2. Por lo tanto: P(Tc k≤tc) = 1 −e−(k 2)tc(2.22) Por estar distribuido exponencialmente, la media y la varianza de estos tiempos de coalescencia Tc kvienen dadas por: E(Tc k) = 2 k(k−1) (2.23) V ar(Tc k) = 2 k(k−1)2 (2.24) De la ecuación (2.23), vemos que según disminuye k, es decir, el número de genes entre los que esperamos que dos encuentren su ancestro común, aumenta el valor esperado de Tc k. Por tanto, el tiempo de coalescencia cuando solamente quedan dos genes para encontrar su ancestro común es el que se espera que sea más largo. Es posible describir un algoritmo que crea genealogías para ngenes, considerando la expresión continua del tiempo. A continuación vemos los distintos pasos que definen el algoritmo. Algoritmo 1 1. Empezar con k=ngenes. 2. Simular el tiempo de espera Tc khasta el siguiente evento de coalescencia, sabiendo que Tc k∼Exp k 2. 24 3. Escoger aleatoriamente una pareja (i, j)con 1≤i<j≤kentre los k 2posibles pares. 4. Convertir iyjen un único gen y disminuir el tamaño de la muestra en una unidad, k→k−1. 5. Si k > 1ir al segundo paso, en otro caso parar. Con los pasos de este algoritmo podemos crear un código en R, como el que se muestra en el Apéndice A.1, para generar un árbol. Así, a continuación presentamos un árbol de una muestra de 5genes, generado con este algoritmo, donde en la izquierda vemos el tiempo continuo, y en la derecha la correspondencia en generaciones. 0.0 1.0 2.0 Algoritmo 1 5 2 4 1 3 4N N 2N 3N Tiempo Figura 2.4: Árbol generado a partir del algoritmo 1. A la izquierda tenemos el tiempo continuo y a la derecha las generaciones a las que se corresponde. 25 2.2.3. Modelo alternativo: Modelo de Moran Ya hemos comentado que la teoría de la coalescencia surge basándose principalmente en el modelo de Wright-Fisher, pero existe otro modelo muy estudiado que también permite la derivación de esta teoría: el modelo de Moran, que fue creado en 1958. La principal razón de su importancia es que, al contrario del modelo de Wright-Fisher, este sí que considera generaciones que se solapan. Además, desde el punto de vista matemático, muchos resultados obtenidos con exactitud bajo el modelo de Moran, eran simplemente aproximaciones en el de Wright-Fisher. Consideremos de nuevo, por simplicidad, una población de 2Nindividuos o genes. La nueva generación se formará a partir de la anterior mediante la selección aleatoria de un gen, para dar lugar a otro gen, y de un gen para morir. El gen que muere no puede ser el gen que va a dar lugar a otro nuevo gen y todos los demás sobreviven a la siguiente generación. Cabe destacar que la forma en la que está construido este modelo impide que más de dos genes encuentren a su ancestro común en la anterior generación. En la Figura 2.5 se percibe mejor la manera en la que se desarrolla el modelo. Los puntos verdes representan genes que mueren y los rosas, genes que dan lugar a otros nuevos. Figura 2.5: Ejemplo del modelo de Moran. La probabilidad de que dos genes encuentren su ancestro común en la anterior generación es 1 N(2N−1). Esto se debe a que tenemos 2N 2=N(2N−1) posibles parejas y solamente una de ellas encuentra su ancestro común. Por lo tanto, el tiempo de espera hasta encontrar el MRCA de dos genes, sigue una distribución geométrica de parámetro 1 N(2N−1) y, entonces, realizando un razonamiento análogo al que se llevó a cabo para el modelo de Wright-Fisher, se puede estudiar este modelo midiendo el tiempo en unidades de N(2N−1) generaciones. Vemos, por tanto, que el análisis de la coalescencia podría ser llevado a cabo también con este modelo, de una manera bastante similar a la empleada en el modelo de WrightFisher. 26 2.3. Medidas del tamaño de una genealogía Volviendo al modelo de Wright-Fisher, que es el que hemos analizado con más precisión, como ya sabemos cómo tiene lugar la coalescencia entre genes de una población, vamos a analizar un par de resultados interesantes relacionados con el tamaño y la forma de una genealogía. El tiempo TMRCA hasta que se encuentra el ancestro común más reciente es de claro interés a nivel biológico a la hora de estudiar una población, así como también lo es otra medida, Ttotal, que es la longitud total de la genealogía. Ttotal también tiene una gran importancia biológica porque indica el tiempo en el que las mutaciones pueden haber ocurrido en la historia de una muestra. Por lo tanto estudiemos desde un punto de vista matemático estas dos cantidades, considerando una muestra de ngenes: el tiempo hasta el ancestro común más reciente de la muestra completa (que al final es simplemente la altura del árbol), TMRCA, y la longitud total de todas las ramas de la genealogía, Ttotal. Como Tc kes el tiempo en la historia de la muestra en el que hubo exactamente klinajes se tiene: TMRCA = n X k=2 Tc k(2.25) Ttotal = n X k=2 kTc k(2.26) En el árbol representado en la Figura 2.4 podemos ver que TMRCA ≈1,98 y que Ttotal = 2Tc 2+ 3Tc 3+ 4Tc 4+ 5Tc 5≈2·0,86 + 3 ·0,28 + 4 ·0,80 + 5 ·0,04 = 5,96. Como TMRCA yTtotal son funciones de variables aleatorias exponenciales, se pueden obtener sus valores esperados simplemente usando las propiedades básicas de la esperanza y recordando que Tc k∼Exp k 2: E[TMRCA] = n X k=2 E(Tc k) = n X k=2 2 k(k−1) = 2 n X k=2 1 k−1−1 k = 2 1−1 2+1 2−1 3+1 3− · · · +1 n−1−1 n = 2 1−1 n−−−→ n→∞ 2(2.27) E[Ttotal] = n X k=2 kE(Tc k) = n X k=2 k2 k(k−1) = 2 n−1 X k=1 1 k(2.28) Cuando el tamaño de la muestra ntiende a ∞, vemos que E[TMRCA]tiende a 2, mientras que E[Ttotal]aumenta sin límite conforme ncrece. 27 34 Capítulo 3 Coalescencia en el modelo de Wright-Fisher con mutaciones 3.1. Modelo de Wright-Fisher con mutaciones En el análisis de datos genéticos reales es esencial tener en cuenta las mutaciones que pueden experimentar los genes. Consideremos entonces, en el modelo de Wright-Fisher explicado en el capítulo anterior, la posibilidad de mutación, es decir, cada gen que se reproduce es susceptible de pasar un proceso de mutación con probabilidad u. Por lo tanto, tenemos que un gen puede ser copiado a su descendencia sin cambios con probabilidad 1−uy, con probabilidad u, puede ser modificado por una mutación. Cabe destacar que según el tipo de datos que estemos considerando uhace referencia a la tasa de mutación por generación, por locus (posición fija en un cromosoma) o por sitio (ubicación en una secuencia de ADN). En la Figura 3.1 vemos tres generaciones del modelo de Wright-Fisher con mutación. En la segunda fila, es decir, en la segunda generación, vemos dos genes que son copias mutadas de genes de la anterior generación. En la generación del presente también se observa un gen que es una copia modificada del de la generación previa. Figura 3.1: Modelo de Wright-Fisher con mutaciones. 35 Siguiendo un determinado linaje vemos que hay una probabilidad ude que la clase del gen padre en la generación tsea distinta de la clase del gen hijo en la generación t+1. Por lo tanto, el número de generaciones que pasan (comenzando desde el presente) hasta que tiene lugar la primera mutación, que denotaremos por TM, sigue una distribución geométrica de parámetro u. Así, la probabilidad de que un linaje experimente su primera mutación t generaciones hacia el pasado es: P(TM=t) = u(1 −u)t−1(3.1) Considerando, como ya se hizo para el modelo de Wright-Fisher sin mutación, que el tiempo es medido en unidades de 2Ngeneraciones, y denotando a este tiempo hasta que se tiene la primera mutación como Tc M, tenemos que, para valores grandes de 2N, podemos aproximar la distribución de Tc Mpor una exponencial. Para ello recordemos que por seguir TMuna distribución geométrica se tiene que P(TM> t) = (1 −u)t. Así: l´ım 2N→∞ P(TM> t) = l´ım 2N→∞ PTM 2N> tc= l´ım 2N→∞ 1−2Nu 2Ntc2N = 1 −e−θtc/2=P(Tc M> tc)(3.2) En esta ecuación hemos definido dos nuevas variables: tc=t 2N, que describe el tiempo en unidades de 2Ngeneraciones, y θ= 4Nu, que es el llamado parámetro de mutación o tasa poblacional de mutación. Este parámetro se puede interpretar como el número esperado de mutaciones que tienen lugar en dos linajes distintos antes de que encuentren su ancestro común. Recordemos que el tiempo esperado para que tenga lugar una coalescencia entre dos genes es 2Ngeneraciones y entonces, en ese tiempo, es de esperar que tengan lugar 2Nu mutaciones en cada linaje, es decir, 4Nu mutaciones en total. Vemos que Tc M sigue una distribución exponencial de parámetro θ 2. A lo largo de este capítulo consideraremos que estamos trabajando con poblaciones de tamaño grande, de forma que podremos considerar la distribución exponencial del tiempo de espera hasta la primera mutación. Por seguir Tc Mesta distribución exponencial de parámetro θ 2, el número Mde mutaciones que tienen lugar en un intervalo de tiempo de duración tcsigue una distribución de Poisson de parámetro θtc 2. La formalización de esta afirmación, relacionada con los procesos de Poisson, se puede encontrar en la página 297 del libro Probability and Measure de Patrick Billingsley [9]. Así, la probabilidad de que tengan lugar mmutaciones en un tiempo tcviene dada por: 36 P(M=m|Tc M=tc) = θtc 2m m!e−θtc 2,con m= 0,1,2, . . . (3.3) Además: E[M|Tc M=tc] = V ar(M|Tc M=tc) = θtc 2(3.4) Debemos tener en cuenta que estas mutaciones que estamos estudiando no aportan ningún tipo de ventaja o desventaja a los linajes, son mutaciones que reciben el nombre de neutras porque no alteran los patrones de reproductividad que hay en una población y, por tanto, son independientes del proceso genealógico. Para estudiar este tipo de mutaciones se han creado varios modelos matemáticos, pero nosotros nos vamos a centrar en uno de ellos: en el modelo de sitios infinitos. 3.1.1. Modelo de sitios infinitos Este modelo apareció en 1969 de la mano del biólogo y matemático japonés Kimura. En él se considera que los genes son simplemente secuencias de ADN. Además, se asume que las mutaciones tienen lugar siempre en una nueva posición (considerando posición como cada uno de los nucleótidos presentes en una secuencia de ADN), es decir, las mutaciones tienen lugar siempre en un nucleótido distinto. Se puede usar este modelo para describir la evolución de las cadenas de ADN muy largas que presentan una baja tasa de mutación en cada posición. Desde el punto de vista biológico esta baja tasa se justifica teniendo en cuenta que, en general, el número de sitios que varían en una muestra de secuencias reales suele ser bastante menor que el número de sitios que son idénticos en todas las secuencias. En el modelo de sitios infinitos siempre habrá uno o dos posibles alelos en una posición de un conjunto de secuencias, pero nunca más, porque cada posición muta como mucho una vez. El modelo también establece que todas las mutaciones que ocurren en algún momento de la historia de la muestra pueden ser recuperadas pues, como solamente puede haber un cambio en un nucleótido específico, si ese cambio tiene lugar se podrá percibir en todo momento. Para entender mejor este modelo se añade la Figura 3.2 que nos muestra el árbol que indica el proceso evolutivo de una secuencia de ADN. Vemos que a lo largo del tiempo tienen lugar cuatro mutaciones (marcadas con puntos negros). Como cada una de ellas ocurre en un nucleótido diferente, vamos marcando con puntos rosas las posiciones diferenciadoras o segregadoras de los nucleótidos donde tienen lugar. Así, tenemos tantas posiciones diferenciadoras como mutaciones hay en la historia de la secuencia. 37 Modelo de sitios infinitos Figura 3.2: Ejemplo del modelo de sitios infinitos. 3.2. Mutaciones y coalescencia Ahora que ya conocemos cómo se pueden introducir las mutaciones en el modelo de Wright-Fisher y que tenemos un modelo para describir el proceso de mutación, como es el modelo de sitios infinitos, vamos a relacionar el proceso de mutación con la coalescencia. Tenemos que el tiempo Tc Mde espera hasta que se tiene una mutación en un linaje dado, sigue una distribución exponencial de parámetro θ 2. Si tenemos klinajes, en cada uno de ellos podremos considerar el tiempo de espera hasta la primera mutación, que en todos los casos sigue una distribución exponencial de parámetro θ 2y son independientes los unos de los otros. En esta situación podemos calcular el tiempo de espera hasta que se tiene una mutación en alguno de los klinajes (Tc Mk), que no será más que el mínimo de los tiempos de espera de cada linaje por separado. Ya hemos visto que dadas dos variables aleatorias XeYque siguen una distribución exponencial de parámetros λyλ0, respectivamente, entonces min(X, Y )∼Exp(λ+λ0). Así, en este caso Tc Mk∼Exp kθ 2. Recordemos además que el tiempo de espera hasta que dos de los klinajes considerados encuentran a su ancestro común, Tc k, sigue una distribución exponencial de parámetro k 2=k(k−1) 2. El tiempo de espera para que ocurra alguno de estos dos eventos (o de coalescencia o 38 de mutación), es una variable aleatoria que es el mínimo de dos variables exponenciales independientes de parámetros kθ 2yk(k−1) 2. Así este tiempo de espera sigue una distribución exponencial de parámetro kθ 2+k(k−1) 2=k(k−1+θ) 2 Ahora que ya sabemos cómo es la distribución del tiempo de espera hasta un evento de coalescencia, hasta un evento de mutación y hasta un evento de cualquiera de los dos tipos, podemos calcular la probabilidad de que el primer evento que tenga lugar en una muestra con klinajes sea de coalescencia o de que sea de mutación. En general, podemos calcular la probabilidad de que un primer evento sea de un determinado tipo considerando que tenemos dos variables exponenciales de parámetros λ1yλ2. Busquemos entonces la probabilidad de que el primer evento sea un proceso de parámetro λ1, es decir, la probabilidad de que el tiempo T1a un evento exponencial (λ1) sea menor que el tiempo T2a un evento exponencial (λ2). Ya hemos visto en la ecuación (1.51) que dadas dos variables exponenciales XeY de parámetros λyλ0respectivamente, se tiene P(X < Y ) = λ λ+λ0. Por lo tanto en esta situación tenemos: P(T1< T2) = λ1 λ1+λ2 (3.5) Vemos que la probabilidad de que el primer evento sea de un tipo determinado, viene dada por la tasa relativa de ese evento. Aplicando esto a nuestro caso, en el que estamos considerando procesos de coalescencia o de mutación, tenemos que la probabilidad de que un evento de coalescencia sea el primero en ocurrir viene dado por: P(coalescencia |coalescencia o mutación) = k(k−1) 2 kθ 2+k(k−1) 2 =k−1 θ+k−1(3.6) mientras que la probabilidad de que una mutación sea el primer evento que ocurra es: P(mutación |coalescencia o mutación) = θ θ+k−1(3.7) Calculemos ahora la distribución del número de eventos que tienen lugar hasta el primer evento de coalescencia (incluyéndolo) entre klinajes. Escribámoslo en forma general, considerando que estamos ante dos procesos descritos por variables exponenciales con parámetros λ1yλ2y sabiendo que la probabilidad de que el siguiente evento sea de “tipo 1” es λ1 λ1+λ2y la de que sea de “tipo 2” es λ2 λ1+λ2. Como los tiempos de espera vienen descritos por variables exponenciales y esta distribución ya hemos visto que carece de memoria, entonces la probabilidad de que suceda un evento de 39 un determinado tipo, una vez ha ocurrido otro, viene dada por las mismas fórmulas que las obtenidas al estudiar cuál es el primer evento. Por lo tanto, los eventos forman una serie de intentos de Bernoulli con probabilidad de éxito λ1 λ1+λ2pues estamos interesados en encontrar el primer evento de tipo 1. Por lo tanto, el número de eventos Bque han tenido lugar cuando el primer evento de tipo 1 ocurre, sigue una distribución geométrica: P(B=b) = λ2 λ1+λ2b−1λ1 λ1+λ2 (3.8) Por lo tanto, la distribución del número de eventos hasta la primera coalescencia (incluyéndola) entre klinajes sigue una distribución geométrica de parámetro k−1 θ+k−1. Es decir, tenemos: P(B=b) = θ θ+k−1b−1k−1 θ+k−1(3.9) 3.3. Algoritmos de generación de genealogías Como ya conocemos las características de los eventos de mutación y de coalescencia, vamos a generar un par de algoritmos que nos permitan simular genealogías en las que tengan lugar estos dos procesos. Comencemos con el que denotaremos por algoritmo 2. Algoritmo 2 1. Empezar con k=ngenes, siendo nel tamaño de la muestra. 2. Considerar una variable exponencial con parámetro k(k−1+θ) 2que determinará cuándo ocurre un evento. 3. Con probabilidad k−1 k−1+θel evento es un evento de coalescencia y con probabilidad θ k−1+θes un evento de mutación. 4. Si ocurre un evento de coalescencia, escoger aleatoriamente un par de genes para encontrar el ancestro común. Actualizar k:k→k−1. 5. Si ocurre un evento de mutación, escoger un linaje para mutar. Dejar el número k invariante. 6. Continuar hasta que kes 1. Este algoritmo es una extensión del algoritmo 1 que ya hemos visto. Para determinar las características de los genes que hay en el presente, se parte del ancestro común y se van viendo las mutaciones que experimenta hasta llegar a cada uno de los genes que hay en la actualidad. Como estamos considerando el modelo de sitios infinitos, cada vez que 40 hay una mutación, el gen, que recordemos que se interpreta como una secuencia de ADN, experimenta una modificación en uno de sus nucleótidos. A continuación, en la Figura 3.3, mostramos un ejemplo de una genealogía de una muestra de 4genes generada con este algoritmo con ayuda de R (ver código en el Apéndice A.2). En ella vemos un evento de mutación y tres eventos de coalescencia. En la derecha vemos la probabilidad que había de que ocurriese un evento del tipo indicado, Algoritmo 2 Tiempo 3 4 1 2 0.0 0.2 0.4 0.6 0.8 1.0 Coalescencia 3 θ + 3 Coalescencia 2 θ + 2 Mutación θ θ + 1 Coalescencia 1 θ + 1 Figura 3.3: Árbol genealógico diseñado con el algoritmo 2. En la derecha se muestra el tipo de evento y la probabilidad de que ocurra en cada caso. Pasemos ahora a otro algoritmo, el algoritmo 3, que genera genealogías similares a las que acabamos de ver pero el procedimiento seguido es ligeramente distinto. Algoritmo 3 1. Simular la genealogía de ngenes de acuerdo con el proceso de coalescencia con parámetro k 2, siendo kel número de linajes considerados en cada etapa. Esta simulación se puede llevar a cabo con el Algoritmo 1. 41 2. Para cada rama considerar el número de mutaciones, M, que viene dado por una distribución de Poisson de parámetro tcθ 2, con tcla longitud de la rama. 3. Para cada rama los tiempos a los que ocurren los Meventos de mutación son escogidos aleatoriamente. Este algoritmo se basa en el hecho de que el número de mutaciones tiene distribución de Poisson. Considera además que las mutaciones pueden ser introducidas en la genealogía una vez que esta ya está generada. Esto se debe al hecho de que las mutaciones que estamos considerando son neutras y no afectan al patrón reproductivo. A continuación, en la Figura 3.4, presentamos un ejemplo de este algoritmo generado con el código creado en R que se puede ver en el Apéndice A.3. Además, resaltamos la idea de que inicialmente se genera la genealogía sin considerar los procesos de mutación y después simplemente se van añadiendo las mutaciones a cada una de las ramas considerando su distribución de Poisson asociada. 2 1 3 4 Algoritmo 3 2 1 3 4 Figura 3.4: Árbol genealógico diseñado con el algoritmo 3. En la izquierda presentamos la genealogía inicial, sin considerar el proceso de mutación y en la derecha tenemos ya introducidas las mutaciones. 42 3.4. Medidas de polimorfismos en una secuencia de ADN Ahora que ya sabemos cómo se pueden modelizar las mutaciones en una genealogía, vamos a emplear el modelo de sitios infinitos para estudiar un concepto esencial en genética: los polimorfismos en una secuencia de ADN. Antes de nada debemos definir qué es un polimorfismo. Los polimorfismos en el ADN son las diferentes secuencias de ADN entre individuos, grupos o poblaciones. Incluyen distintos niveles de variación, desde un único cambio en la base nitrogenada de un nucleótido, cambios en varias bases o cambios en secuencias. Aquí vamos a considerar los polimorfismos de un único nucleótido (SNP), es decir, aquellos en los que se produce una única variación en la base del nucleótido. [10] Una de las formas más simples de estudiar los polimorfismos del ADN es obteniendo el número de sitios segregadores de una muestra, que denotaremos por S. Los sitios segregadores no son más que las posiciones de los nucleótidos en los que se produce la mutación. Para comprender de formas más visual lo que son los sitios segregadores, presentamos en la Figura 3.5 cinco secuencias en las que se observan en color tres sitios segregadores, es decir, S= 3. Secuencia 1: A CG C TAGTCA Secuencia 2: A GG C TAGTCA Secuencia 3: A GG C TAGTCT Secuencia 4: A GG C AAGTCT Secuencia 5: A CG C AAGTCT Figura 3.5: Secuencias de ADN con los sitios segregadores señalados. Recordemos que en el modelo de sitios infinitos, cada mutación ocurre en una única posición, por lo que toda mutación que ocurra en la evolución de la muestra, será un sitio segregador. Así, el número Sde sitios segregadores en una muestra de tamaño nes igual al número de mutaciones en la evolución de la muestra en el modelo de sitios infinitos. Tenemos que tener en cuenta que estamos considerando una genealogía de longitud total Ttotal. Entonces, el número de mutaciones en una genealogía de esta longitud sigue una distribución de Poisson de parámetro θTtotal 2. Como conocemos la función de densidad de Ttotal: fTtotal (tc) = n X k=2 (−1)kn−1 k−1k−1 2e−k−1 2tc(3.10) 43 Tras esta publicación, en 1998, Nordborg [14] decidió emplear los fundamentos de la teoría de la coalescencia para estudiar la relación que tienen los seres humanos actuales y los neandertales. La muestra neandertal se considera que fue datada en un tiempo ts que oscila entre los 30000 y los 100000 años. Como el ADN mitocondrial se transmite de madres a hijos, asumiendo que la población de mujeres tenía un tamaño de 3400 y que cada generación tiene una vida de unos 20 años, entonces tstoma valores entre 0,44 y1,47 en la escala de tiempo de coalescencia. Debemos mencionar que en todo este apartado vamos a considerar que el tiempo es continuo y trabajaremos con el formalismo continuo de la coalescencia. Lo que realmente quería analizar Nordborg con su estudio era la hipótesis nula de que, cuando los seres humanos y los neandertales coexistían, había emparejamiento aleatorio entre ellos. Dos de los factores que Krings recalcó y que hicieron a Nordborg creer que esta hipótesis era falsa eran, por un lado que Trfuese 4veces mayor que Te, y por otro lado, que la estructura del árbol fuese tal que el neandertal no se juntase con ningún linaje humano hasta que solamente quedaba uno. En un principio Nordborg se planteó la idea de que esta forma del árbol ya permitiese rechazar la hipótesis nula, pero observó que no era una condición suficiente pues si el neandertal y el último linaje humano encuentran su ancestro común en un tiempo muy pequeño no podríamos concluir que no hubiese emparejamiento aleatorio entre ellos. Fue por esto por lo que se planteó considerar también el hecho de que Trfuese mayor o igual que 4Te. Así, con todo esto, Nordborg consideró como p-valor para analizar esta hipótesis nula, la probabilidad de que se diese un modelo tan extremo o más que el establecido por Krings, es decir, la probabilidad de que se diese un árbol con la forma establecida anteriormente y de que Tr≥4Te, lo que escribiremos como: P(árbol y Tr≥4Te)(4.1) En la Figura 4.1 vemos representado el árbol que hace referencia a la situación que estamos considerando, en donde tenemos 986 muestras de humanos actuales y una de un neandertal, que aparece en un tiempo ts. Tenemos también representado el tiempo de coalescencia de todos los humanos (Te) y el de los humanos y el neandertal (Tr). 50 ts Te Tr 986 humanos actuales Neandertal Figura 4.1: Árbol que representa la situación entre las muestras de los humanos y la muestra neandertal. Usaremos An(t)para denotar el número de linajes ancestrales que existen en el tiempo t en el pasado, de una muestra de tamaño ntomada en el presente (en nuestro caso n= 986). Nordborg razonó que, conocido el valor de An(t), la probabilidad de que el árbol tenga esa forma y de que Tr≥4Teson independientes. Por lo tanto podemos escribir: P(árbol y Tr≥4Te) = 986 X k=1 P(árbol|An(ts) = k)P(Tr≥4Te|An(ts) = k)P(An(ts) = k)(4.2) Lo primero que haremos será calcular la probabilidad de que en un tiempo t, partiendo de una muestra de tamaño n, se tengan klinajes, es decir, trataremos de calcular gn,k(t) = P(An(t) = k). Si nos fijamos, esta probabilidad es la misma que la de que antes del tiempo t, hayan tenido lugar exactamente n−keventos de coalescencia. Para calcular gn,k haremos un estudio ligeramente distinto si nos encontramos en el caso k= 1 de los demás. Comencemos con gn,1(t). En este caso estamos buscando la probabilidad de que hayan tenido lugar n−1eventos de coalescencia antes del tiempo t, es decir, la probabilidad de que todos los individuos de la muestra de tamaño nhayan encontrado a su ancestro común más reciente. Sabiendo esto y recordando que ya hemos obtenido la distribución de TMRCA 51 (ver ecuación (2.41)), se tiene: gn,1(t) = Zt 0 fTMRCA (x)dx =Zt 0 n X i=2 i 2e−(i 2)x n Y j=2,j6=ij 2 j 2−i 2dx = n X i=2 1−e−(i 2)tn Y j=2,j6=ij 2 j 2−i 2(4.3) En este procedimiento simplemente hemos sustituido el valor de fTMRCA , que ya habíamos calculado en apartados anteriores, y hemos empleado la integral de una función exponencial. Pasemos ahora al caso 2≤k≤(n−1). En esta situación, para calcular gn,k(t)es necesario tener en cuenta que el (n−k)-ésimo evento de coalescencia ocurre antes de t, pero, sin embargo, el (n−k+ 1)-ésimo ocurre después de t. Vamos a considerar una nueva variable: Tn,k =Pn i=k+1 Tc i, que denota el tiempo que pasa hasta que tiene lugar el (n−k)-ésimo evento de coalescencia. Una interpretación equivalente de esta variable es que describe el tiempo que pasa desde que tenemos nlinajes hasta que tenemos k. La distribución de esta nueva variable se obtiene de la misma manera que obtuvimos fTMRCA en su momento, considerando la convolución de variables aleatorias exponenciales independientes. Así, llegamos, para 2≤k < (n−1), a que: fTn,k (t) = n X i=k+1 i 2e−(i 2)t n Y j=k+1,j6=ij 2 j 2−i 2(4.4) En el caso k=n−1tenemos que Tn,n−1no es más que el tiempo en el que hay n linajes, es decir, no es más que Tc n, y por tanto fTn,n−1(t) = n 2e−(n 2)t. Para calcular gn,k(t)simplemente tendremos que considerar que el (n−k)-ésimo evento de coalescencia ocurre en un instante xantes de ty que el tiempo en el que hay klinajes (Tk) es mayor que t−xpues el (n−k+ 1)−ésimo evento de coalescencia tiene que ocurrir después de t. Esta idea se plasma en la flecha temporal de la derecha. xMomento en el que ocurre el n−k evento de coalescencia t Tk En esta situación, para 2≤k≤(n−1) tenemos que: gn,k(t) = Zt 0 fTn,k (x)Z∞ t−x fTc k(y)dydx (4.5) 52 Así, en el caso 2≤k < (n−1) simplemente sustituyendo las funciones de densidad de las variables, integrando y haciendo una serie de cálculos llegamos a la siguiente expresión: gn,k(t) = n X i=ki 2 k 2e−(i 2)t n Y j=k,j6=ij 2 j 2−i 2(4.6) Para k=n−1, integrando llegamos a: gn,n−1(t) = n 2 n−1 2−n 2he[(n−1 2)−(n 2)]t−1ie−(n−1 2)t(4.7) Para el caso k=n,gn,k(t)indica la probabilidad de que cuando estemos en el tiempo taún tengamos los nlinajes iniciales, es decir, que no se haya producido ningún evento de coalescencia en todo ese tiempo. Esto es sinónimo a que el tiempo Tc nsea mayor que t. Así: gn,n(t) = Z∞ t fTc n(x)dx =e−(n 2)t(4.8) De esta forma, ya conocemos la probabilidad de que en el tiempo thaya klinajes, para cualquier valor de k, partiendo de nlinajes iniciales, es decir, ya conocemos gn,k(t). Pasemos ahora a calcular la probabilidad relacionada con la forma del árbol. Esta probabilidad hace referencia al hecho de que los neandertales no comparten ningún ancestro común con el ser humano antes de que haya un único linaje de los humanos. Así, calcularemos esta probabilidad condicionada al hecho de que A986(ts) = k. Por lo tanto, buscamos la probabilidad de que un linaje particular (el linaje del neandertal), en una muestra de tamaño k+ 1 (considerando los klinajes humanos presentes en el tiempo tsy el linaje neandertal) no encuentre el ancestro común con ningún otro linaje hasta el final, cuando solamente queden dos linajes. Sabemos que en general si tenemos jlinajes, hay j 2=j(j−1) 2 posibles parejas que pueden encontrar su ancestro común. Dentro de ellas j−1involucrarán al linaje del neandertal. Así, hay j 2−(j−1) posibles parejas de linajes que pueden encontrar su ancestro común de entre las j 2parejas existentes. Entonces, la probabilidad que buscamos será simplemente el producto desde j= 3 (cuando hay 2linajes humanos y uno neandertal) hasta j=k+ 1, (cuando hay klinajes humanos y uno neandertal) de los eventos de coalescencia que no involucran al linaje del neandertal entre los posibles totales que se podrían dar en general: (j 2)−(j−1) (j 2). P(árbol|A986(ts) = k) = k+1 Y j=3 j 2−(j−1) j 2!= k+1 Y j=3 1−(j−1) j(j−1) 2! = k+1 Y j=3 1−2 j= k+1 Y j=3 j−2 j =1·2· · · (k−3)(k−2)(k−1) 3·4· · · (k−2)(k−1)k(k+ 1) =2 k(k+ 1) (4.9) 53 Por último nos queda calcular la probabilidad de que Trsea mayor o igual que 4veces Te, sabiendo que en el tiempo tsse tienen klinajes humanos. Al igual que hicimos antes para calcular gn,k vamos a distinguir el caso k= 1 de los demás. Comencemos para k≥2. Así: P(Tr≥4Te|An(ts) = k) = P(Tr−4Te≥0|An(ts) = k) =P((Tr−ts)−4Te≥ −ts|An(ts) = k) =P(Tk+1,1−4Tn,1≥ −ts) =PTk+1,1 4−Tn,1≥−ts 4(4.10) En este razonamiento hemos restado tsa ambos lados de la desigualdad y hemos usado que Tr−ts=Tk+1,1, pues es el tiempo que pasa entre que tenemos k+ 1 linajes (el neandertal y los khumanos) hasta que solamente tenemos 1. Además también hemos empleado que Te=Tn,1, ya que esta variable determina el tiempo en el que se pasa de los nlinajes humanos iniciales a un único linaje humano. Finalmente hemos dividido todo entre 4. En la Figura 4.2 tenemos el árbol con k= 2 y vemos marcado Tr−ts. ts Te Tr 986 humanos actuales Neandertal Tr−ts Figura 4.2: Árbol que representa la relación entre humanos y el neandertal junto con Tr−ts señalado. Para obtener la distribución de Tk+1,1 4calculamos su función de distribución y derivamos con el fin de obtener la función de densidad. Así llegamos a: fTk+1,1 4 (t) = k+1 X i=2 4i 2e−4(i 2)t k+1 Y j=2,j6=ij 2 j 2−i 2(4.11) 54 La distribución de Tn,1ya la hemos visto antes en su fórmula general: fTn,1(t) = n X i=2 i 2e−(i 2)t n Y j=2,j6=ij 2 j 2−i 2(4.12) Las dos variables TeyTk+1,1son independientes, pues estamos considerando por un lado a los humanos y por otro a los humanos más el neandertal. Así, ahora que ya sabemos todo esto, podemos calcular el valor de la probabilidad buscada. Para ello consideramos la siguiente expresión (en la que ya hemos tenido en cuenta que las dos variables son independientes) y empleamos el cambio de variable y=z−x: P(Tr≥4Te|An(ts) = k) = Z∞ ts 4"Z∞ x−ts 4 fTn,1(x)fTk+1,1 4 (z)dz#dx +Zts 4 0Z∞ 0 fTn,1(x)fTk+1,1 4 (z)dzdx = n X i=2 k+1 X r=2 "1−4r 2 i 2+ 4r 2e−(i 2)ts 4#pipr(4.13) donde: pi= n Y j=2,j6=ij 2 j 2−i 2pr= k+1 Y s=2,s6=rs 2 s 2−r 2(4.14) Para el caso k= 1 se tiene que Tk+1,1=T2,1=Tc 2. Como la distribución de Tc 2sabemos que es exponencial de parámetro 1, podemos calcular fácilmente la probabilidad buscada. Así: P(Tr≥4Te|k= 1) = Z∞ ts 4"Z∞ x−ts 4 fTn,1(x)fT2 4 (z)dz#dx = n X i=2 "1−4 i 2+ 4e−(i 2)ts 4#n Y j=2,j6=ij 2 j 2−i 2(4.15) Ahora ya tenemos todas las componentes necesarias para calcular el p-valor deseado. De esta forma, sustituyendo con ayuda de R (ver código en el Apéndice A.5) para n= 986 y recordando que consideramos que la población de mujeres era de 3400 y que cada generación duraba 20 años, podemos calcular las siguientes cantidades: ts(en años) 30000 100000 E[A986(ts)] 4,87 1,75 P(árbol)0,085 0,56 P(árbol y Tr≥4Te)0,0020 0,023 55 Vemos que E[A986(ts)] toma valores bastante pequeños para los dos valores de tsconsiderados. Esto indica que el número esperado de linajes humanos que siguen presentes en un tiempo tses bastante pequeño (comparado con los 986 linajes que había inicialmente). Además observamos que el p-valor que hemos definido toma valores suficientemente pequeños (menores que 0,05) para los dos tsconsiderados. Por lo tanto, teniendo esto en cuenta, podemos rechazar la posibilidad de que haya habido emparejamiento aleatorio entre los seres humanos y los neandertales al nivel del 0,2 % (en el caso de ts= 30000) y al nivel del 2,3 % (en el caso de ts= 100000). Desde el año 1998 hasta la actualidad el estudio de los antepasados del ser humano moderno ha avanzado significativamente. Hoy en día se tiene mucha más información sobre este tema y han sido encontradas más muestras paleontológicas de homínidos antiguos. Además, se ha conseguido secuenciar al completo el genoma humano y el ADN de varios individuos neandertales. Con todo esto, la mayoría de los expertos en el tema establecen que sí que hay contribución de los neandertales en los seres humanos actuales ([15]), mientras que sigue habiendo alguno ([16]) que basándose en la comparación de ADN mitocondrial, establece que en el caso de que hubiese contribución esta sería pequeña. En general, la tendencia hoy en día, con todos los datos que se tienen, es la de creer que sí que hubo un momento en el que los neandertales y los humanos convivieron y se emparejaron entre sí. 56 Apéndice A Códigos de R A.1. Algoritmo 1 ##Algoritmo 1 #Número de genes n=5 # Creación del cluster dd <- dist(scale(seq(1:n)), method ="euclidean") hc <- hclust(dd, method ="ward.D2") hc$labels=c(1:n) hc$order=c(1:n) #Inicialización de variables x=c(-n:-1) i=1 T_k=0 #Algoritmo while(length(x)>1){ lambda=choose(length(x),2) T_k=rexp(1,rate=lambda)+T_k hc$height[i]=T_k c=sample(x,2,replace=F) 57 x=setdiff(x,c) x=append(x,i,after=0) hc$merge[i,1]=c[1] hc$merge[i,2]=c[2] i=i+1 } #Representación del árbol hcd <- as.dendrogram(hc) plot(hcd, type ="rectangle",main="Algoritmo 1") A.2. Algoritmo 2 ##Algoritmo 2 #Número de genes n=6 # Creación del cluster dd <- dist(scale(seq(1:n)), method ="euclidean") hc <- hclust(dd, method ="ward.D2") hc$labels=c(1:n) hc$order=c(1:n) #Inicialización de variables x=c(-n:-1) i=1 T_k=0 theta=2 k=n tm=c() #Algoritmo while(length(x)>1){ 58 lambda=k*(k-1+theta)/2 T_k=rexp(1,rate=lambda)+T_k r=runif(1,0,1) l1=(k-1)/(k-1+theta) if(r<=l1){ hc$height[i]=T_k c=sample(x,2,replace=F) x=setdiff(x,c) x=append(x,i,after=0) hc$merge[i,1]=c[1] hc$merge[i,2]=c[2] i=i+1 k=k-1 } else{ tm=c(tm,T_k) } } # Gráfico hcd <- as.dendrogram(hc) plot(hcd, type ="rectangle",ylab ="Tiempo",main="Algoritmo 2") #Dibujar los puntos de las mutaciones library(ggdendro) segmentos=dendro_data(hcd)$segments d=c() if (length(tm)>0){ for (i in 1:length(segmentos$x)){ if (segmentos$y[i]==segmentos$yend[i]){ d=append(d,i) } } segmentos=segmentos[-d,] m=nrow(segmentos) n1=c() 59 66 Bibliografía [1] Sheldon Ross. A First Course in Probability. Pearson, 2014. [2] Norman Lloyd Johnson, Samuel Kotz, and Narayanaswamy Balakrishnan. Discrete multivariate distributions, volume 165. Wiley New York, 1997. [3] Anirban DasGupta. Fundamentals of probability: a first course. Springer Science & Business Media, 2010. [4] Marco Taboga. Joint moment generating function.https://www.statlect.com/ fundamentals-of-probability/joint-moment-generating-function, 2017. [5] Jotun Hein, Mikkel Schierup, and Carsten Wiuf. Gene Genealogies, Variation and Evolution: A Primer in Coalescent Theory. Oxford University Press, USA, 2004. [6] John Wakeley. Coalescent Theory: An Introduction. Roberts and Company Publishers, 2009. [7] José Manuel Sánchez Muñoz. El problema de Basilea. Lecturas matemáticas, 35(2):199–228, 2014. [8] Simon Tavaré, David J Balding, Robert C Griffiths, and Peter Donnelly. Inferring coalescence times from DNA sequence data. Genetics, 145(2):505–518, 1997. [9] Patrick Billingsley. Probability and measure. John Wiley & Sons, 1995. [10] Yamin Liu. Genetic Diversity and Disease Susceptibility. BoD–Books on Demand, 2018. [11] GA Watterson. On the number of segregating sites in genetical models without recombination. Theoretical population biology, 7(2):256–276, 1975. [12] G Danesh, B Elie, and S Alizon. Early phylodynamics analysis of the covid-19 epidemics in france using 194 genomes. Technical report, 2020. 67 [13] Matthias Krings, Anne Stone, Ralf W Schmitz, Heike Krainitzki, Mark Stoneking, and Svante Pääbo. Neandertal DNA sequences and the origin of modern humans. cell, 90(1):19–30, 1997. [14] Magnus Nordborg. On the probability of neanderthal ancestry. American journal of human genetics, 63(4):1237, 1998. [15] Anders Bergström, Chris Stringer, Mateja Hajdinjak, Eleanor ML Scerri, and Pontus Skoglund. Origins of modern human ancestry. Nature, 590(7845):229–237, 2021. [16] David Serre, André Langaney, Mario Chech, Maria Teschler-Nicola, Maja Paunovic, Philippe Mennecier, Michael Hofreiter, Göran Possnert, and Svante Pääbo. No evidence of neandertal mtdna contribution to early modern humans. In Early modern humans at the Moravian gate, pages 491–503. Springer, 2006. 68