El problema

Entre els problemes que preteníem resoldre en classe, ens ha sorprés el següent:

“De les 100 diferències en llocs de restricció entre dues soques de l’escarabat de la farina, Tribolium castaneum, que es creuen i després es deixen aparellar a l’atzar, quin nombre de llocs de restricció esperaries que continuaren segregant després de 10 generacions assumint una mida poblacional efectiva de 80 individus? Quants esperaries que no s’hagueren fixat després de 50 generacions?”

És un cas de canvi de freqüències al·lèliques com a conseqüència de la deriva genètica. Els llocs de restricció actuen com a marcadors moleculars: són loci amb un al·lel diferent en cada llinatge, només un dels quals és reconegut i tallat per un enzim de restricció. En la generació híbrida tots els individus són heterozigots i les freqüències al·lèliques són 0.5 en tots els loci.

El fet que tots els individus siguen heterozigots és conseqüència del creuament experimental, i implica un coeficient d’endogàmia negatiu. En tan sols una generació d’aparellament aleatori, les freqüències genotípiques s’aproximaran molt a les esperades en l’equilibri de Hardy-Weinberg, i per tant $H_0 = 2pq = 0.5$. Si la població fóra infinita, l’heterozigositat es mantindria a eixe nivell exactament.

La deriva

La deriva genètica és un procés estocàstic. És impossible predir la trajectòria concreta que seguirà la freqüència d’una al·lel en una població concreta. Sabem, però, que a la llarga un dels dos al·lels presents en un locus desapareixerà i l’altre es fixarà. Si repetírem el procés amb un nombre elevat de poblacions, a vegades es fixaria un al·lel, a vegades, l’altre. Ocasionalment (per exemple, quan acaba fixant-se l’al·lel menys freqüent) l’heterozigositat pot augmentar, encara que només per tornar a disminuir després. L’equació $H_t = H_0\left( 1 - \frac{1}{2N_e}\right) ^t$ descriu la disminució de l’heterozigositat mitjana entre totes estes trajectòries possibles. No representa necessàriament la trajectòria de l’heterozigositat en cap població concreta.

Així doncs, en qualsevol moment existeix la possibilitat de què un al·lel desaparega i l’heterozigositat, per tant, esdevinga zero en una població com a conseqüència de la deriva: fins i tot a pesar de què $H_t > 0$. El que ens interessa saber és la probabilitat de què es fixe un al·lel (o l’altre) en cada generació que passa. Conéixer la mitjana d’una distribució ($H_t$) no és suficient per saber quina és la probabilitat que la distribució assigna al valor 0.

L’endogàmia

De passada, convé aclarir una possible confusió: si l’aparellament és aleatori, i per tant les freqüències genotípiques estan en equilibri de Hardy-Weinberg, per què es diu que la deriva fa augmentar el coeficient d’endogàmia $F$? De fet, la fòrmula $(1 - F_t)=(1-F_0)\left( 1 - \frac{1}{2N_e}\right) ^t$ indica un augment de l’endogàmia a cada generació que passa. I també és cert que en tot moment una població concreta amb aparellament aleatori no té cap excés d’homozigots respecte del que cabria esperar amb les freqüències al·lèliques d’eixe moment.

La solució d’aquesta aparent contradicció és que les freqüències al·lèliques de referència no són les actuals en cada moment, sinó les originals. En tota població sotmesa a deriva s’acumula un excése d’homozigots respecte dels que hi hauria si les freqüències al·lèliques s’haguessin quedat com al principi.

Una altra manera de veure-ho és la següent. El coeficient d’endogàmia, F, mesura la probabilitat de què en un zigot es reunisquen dos al·lèls idèntics per descendència, és a dir que són còpies d’un mateix avantpassat comú relativament recent. Com de recent ha de ser l’avantpassat comú de dues còpies per considerar-les idèntiques per descendència és una qüestió pragmàtica, i normalment es tria el moment de fundació d’una població, el moment de subdivisió, etc. És a dir, el moment a partir del qual actua més la deriva. En el cas del problema, és el moment en què la descendència del primer creuament comença a reproduir-se i les freqüències al·lèliques eren 0.5. En qualsevol altre moment posterior, les freqüències al·lèliques hauran canviat, i adoptar les noves freqüències com a referència seria equivalent a considerar idèntiques per descendència únicament les còpies descendents del mateix ancestre en la generació immediatament anterior, cosa que no tindria cap significat pràctic.

La dificultat

Predir l’heterozigositat mitjana al llarg d’un procés de deriva és fàcil. Però predir la probabilitat de cada freqüència al·lèlica possible en cada una de les generacions al llarg del procés és molt més difícil. Aquest és un problema clàssic al qual Fisher, Haldane, Wright i altres donaren solucions parcials en la primera meitat del segle XX. No va ser fins el 1955 que es va trobar una solució analítica satisfactòria. Kimura (Kimura 1955) va haver d’introduir les equacions de difusió a la genètica de poblacions per tal de resoldre-ho.

En una població diploide de mida $N$, les freqüències al·lèliques només poden prendre $2N + 1$ valors diferents, entre 0 i 1, corresponents al nombre de còpies de l’al·lel que poden haver en la població. Tanmateix, l’aproximació feta amb equacions de difusió representa la freqüència al·lèlica com una variable contínua.

Les equacions de difusió estan fora dels objectius del programa, lamentablement. Arribats a este punt, pot semblar un fracàs no poder respondre les preguntes de l’enunciat del problema. Al contrari: reconéixer la complexitat d’un problema i les limitacions pròpies és necessari per poder donar, algun dia, la resposta correcta.

Les cadenes de Markov

Imagine que a Kimura li hagués encantat disposar dels ordinadors d’avuí en dia i de llenguatges de programació tan pràctics com R. Amb aquestes eines podem trobar solucions numèriques al problema de les probabilitats de cada freqüència al·lèlica possible en cada generació, de forma molt més senzilla que en aquells temps.

La deriva és un procés de Markov, perquè la probabilitat de què la freqüència d’un al·lel passe de $i$ (en número de còpies) a $j$ en una generació no depén de com s’han arribat a tenir $i$ còpies en la generació anterior. És un procés sense memòria, i es pot representar mitjançant una matriu de transicions $\mathbf{P}$ que especifica la probabilitat de passar a tenir $j$ còpies d’un al·lel en una generació, sabent que en l’anterior n’hi havia $i$. Anomenant $\mathbf{X}_t$ al vector de probabilitats dels valors possibles de freqüència al·lèlica en el moment $t$:

$$\begin{align} \mathbf{X}_{t+1} &= \mathbf{X}_t * \mathbf{P}\\ \mathbf{X}_t &= \mathbf{X}_0 * \mathbf{P}^t \end{align}$$

Els vectors $\mathbf{X}_t$ i $\mathbf{X}_{t+1}$, indiquen les probabilitats de què l’al·lel en qüestió tinga freqüència 0, $\frac{1}{2N}$, $\frac{2}{2N} \ldots 1$. En la matriu $\textbf{P}$, l’element de la fila $i$ i la columna $j$ representa la probabilitat de passar de $i$ a $j$ còpies d’aquell al·lel en una generació, la qual ve determinada per la distribució binomial: $P_{ij} = {2N\choose j} \left( \frac{i}{2N}\right) ^j\left( 1 - \frac{i}{2N}\right) ^{2N-j}$.

En elevar la matriu $\mathbf{P}$ a $t$, obtenim una nova matriu de transició que indica les probabilitats de passar de $i$ a $j$ no en una sinó en $t$ generacions.

En realitat tenim tots els números necessaris per fer aquests càlculs, i en R es pot implementar de la manera següent

library(expm)                 # Llibreria necessària per elevar una matriu a una potència.
library(RColorBrewer)
library(plyr)

N <- 80                       # Mida poblacional efectiva.
X <- rep(0, 2*N + 1)          # Hi ha 2N+1 valors possibles de freqüència al·lèlica: 0, 1/2N, 2/2N... 1.
                              # Els iniciem tots a zero, i donem tota la probabilitat al valor en posició
X[N + 1] <- 1                 # N + 1: hi ha N còpies de l'al·lel, que fan una freqüència de 0.5.

P <- matrix(                  # Creem una matriu i l'omplim, fila a fila, amb els valors de la
  dbinom(                     # probabilitat bionmial corresponents a:
    rep(0:(2 * N), 2 * N + 1),     # obtenir "i" còpies de l'al·lel, entre 0 i 2N, en 2N+1 columnes,
    2 * N,                         # d'un total de 2N còpies, on...
    rep(0:(2 * N), each = 2 * N + 1) / (2 * N)),  # ...la freqüència de l'al·lel en la generació anterior és j/2N,
                                                  # repetint cada valor per totes les columnes d'una fila.
  ncol = 2 * N + 1,           # El nombre de columnes de la matriu (les files s'ajusten adientment).
  byrow = TRUE)               # La matriu s'ompli per files.

# Moments del temps, en generacions, on volem inspeccionar la distribució de freqüències:
Generacions <- c(N/10, N/5, N/2, N, 2*N, 3*N)
Densitats <- sapply(Generacions,
                    function(t) X %*% (P %^% t))

Colors <- brewer.pal(length(Generacions), name='RdYlGn')
plot(c(0,1), c(0, max(Densitats[,1])),
     xlab = 'Freqüència al·lèlica',
     ylab = 'Probabilitat',
     type = 'n'
     )
l_ply(1:length(Generacions), function(x) {
  lines((0:(2 * N)) / (2 * N), Densitats[,x], col = Colors[x])
})

El color de les línies indica el nombre de generacions que han passat, des del més roig (8) fins el més verd (240). Aquestes formes reprodueixen molt bé les que va obtenir Kimura amb les equacions de difusió. A mesura que passa el temps, la distribució s’aplana pel centre i s’eleva en els extrems.

L’heterozigositat

Coneixent les probabilitats de les freqüencies al·lèliques en cada moment, podem calcular l’heterozigositat esperada com a $H_t=\sum_{p=0}^{1}f_t(p)\cdot 2p(1-p)$, on $f_t(p)$ és la probabilitat de què la freqüència al·lèlica siga $p$ en el moment $t$. Aquesta $H_t$, estimada a partir de la recursió de les probabilitats binomials, hauria de coincidir amb la fòrmula que coneixem: $H_t=H_0\left( 1 - \frac{1}{2N}\right) ^t$. Comprovem-ho.

Heterozigositat <- sapply(1:length(Generacions),
                          function(t) {
                            p <- 0:(2 * N) / (2 * N)
                            q <- 1 - p
                            return(sum(Densitats[,t] * 2 * p * q))  # Heterozigositat mitjana esperada en t.
                          })
plot(c(1,160), c(0.5*(1-1/(2*N))^160, 0.5), type='n', xlab='Generacions', ylab='Heterozigositat')
points(Generacions, Heterozigositat)                   # punts calculats
lines(1:160, 0.5 * (1 - 1/(2*N))^(1:160), col='red')   # corba teòrica

També podem representar la distribució de probabilitats de l’heterozigositat al llarg del temps. A cada valor de freqüència al·lèlica, $p$, li correspon un valor de $H=2p(1-p)$, i cada valor d’heterozigositat es pot obtenir de dues maneres, intercanviant $p$ per $1-p$. Per tant, la probabilitat de cada valor de $H$ serà la suma de les probabilitats dels valors de $p$ que li corresponen.

DensiHet <- Densitats[1:N,] + Densitats[(2 * N + 1):(N + 2),]
DensiHet <- rbind(DensiHet, Densitats[N+1,])   # Densitats de p=q=0.5
H <- 2 * (0:N)/(2 * N) * (1 - (0:N)/(2 * N))   # Valors d'heterozigositat de cada fila de DensiHet
plot(H, DensiHet[,1], type='n', xlab='Heterozigositat', ylab='Probabilitat')
l_ply(1:length(Generacions), function(x) {
  lines(H, DensiHet[, x], col=Colors[x])
})

Resposta al problema

Aleshores, quants dels 100 loci inicialment polimòrfics dels quals parla el problema esperem que s’hagen fixat per a un al·lel o per a l’altre al cap de 10 o 50 generacions? En la mesura en que els 100 loci siguen independents (no estiguen lligats), les trajectòries de les seues freqüències al·lèliques seran independents, com si es tractaren de 100 poblacions diferents. Per tant, la fracció d’aquells loci que deixaran de ser polimòrfics és igual a la probabilitat de fixació d’un o altre al·lel: $P(\mathrm{fixació}|t=10) = P(i=0 | t=10) + P(i=1 | t=10)$.

P_fixacio_10 <- (X %*% (P %^% 10))[1] + (X %*% (P %^% 10))[2 * N + 1]
P_fixacio_10

## [1] 5.294322e-08

P_fixacio_50 <- (X %*% (P %^% 50))[1] + (X %*% (P %^% 50))[2 * N + 1]
P_fixacio_50

## [1] 0.03554453

Així doncs, dels 100 loci inicialment polimòrfics, esperem que al cap de 10 generacions d’aparellament aleatori se n’hagen fixat 100 * P_fixacio_10 = 0. I al cap de 50 generacions, 100 * P_fixacio_50 = 3.55.

Equacions de difusió

Un altre dia parlarem d’elles.