Simulering og inferens i matematiske modeller fra biologien
Lars Nørvang Andersen
Institut for Matematik
Aarhus, Marts 2015
Gorilla projekt - data
Motivationen for projektet var at studere artsdannelse i
forbindelse kortlægningen af gorillaens genom.
Model
t=TA t=TA t=TA
t=0 Oprindelig
population
Population 1
Population 2
Fortid
Nutid m1→2
m2→1 m1→2 m2→1
Fokus lå på en beskrivelse af den proces som ledte frem til de to populationer, som observeres. Specifikt:
Tidspunktet for opsplitning af den oprindelige population.
Tilstedeværelse/fravær af migration.
Gorilla projekt - data
Som data betragtes kolonner i DNA alignments.
Gorilla projekt - data
Datasættet kan således opsummeres i en 10 × 10 tabel:
Eastern Gorilla
AA CC GG TT AC AG AT CG CT GT
WesternGorilla
AA nAAAA nAACC nAAGG nAATT nAAAC nAAAG nAAAT nAACG nAACT nAAGT CC nCCAA nCCCCnCCCC nCCGG nCCTT nCCAC nCCAG nCCAT nCCCG nCCCT nCCGT GG nGGAA nGGCC nGGGG nGGTT nGGAC nGGAG nGGAT nGGCG nGGCT nGGGT TT nTTAA nTTCC nTTGG nTTTT nTTAC nTTAG nTTAT nTTCG nTTCT nTTGT AC nACAA nACCC nACGG nACTT nACAC nACAG nACAT nACCG nACCT nACGT AG nAGAA nAGCC nAGGG nAGTT nAGAC nAGAG nAGAT nAGCG nAGCT nAGGT AT nATAA nATCC nATGG nATTT nATAC nATAG nATAT nATCG nATCT nATGT CG nCGAA nCGCC nCGGG nCGTT nCGAC nCGAG nCGAT nCGCG nCGCT nCGGT CT nCTAA nCTCC nCTGG nCTTT nCTAC nCTAG nCTAT nCTCG nCTCT nCTGT GT nGTAA nGTCC nGTGG nGTTT nGTAC nGTAG nGTAT nGTCG nGTCT nGTGT
For at fortsætte, skal vi opstille en statistisk model for data således at vi udlede likelihood-funktionen L(T A , m 1→2 , m 2→1 ).
Vores model skal sætte os i stand til at tilskrive
sandsynligheder til cellerne i vores tabel, hvorefter vores
likelihood-funktionen kan bruges til estimere parametrene.
Wright-Fisher modellen - to alleler
Nutid
Fremtid
Antager blandt andet: Konstant populationsstørrelse, ingen
mutationer, ikke-overlappende generationer.
R
R er et open source statistikprogram, som kan hentes på adressen:
http://www.r-project.org/
Et populært (og open source) udviklingsmiljø er programmet RStudio, som kan findes på adressen:
http://www.rstudio.com/
R
R-kode til at generere een generation:
gen = c('A','A','A','A','a','a','a','a') Liste med 4 A-er og 4 a-er sample(gen,replace = TRUE) sample bruges til udtræk med tilbagelægning
[1] "a" "a" "a" "a" "a" "a" "A" "a"
Vi kan iterere ovenstående og gemme antallet af ’A’-er i hver generation, indtil der kun er ’A’-er eller ’a’-er.
gen=c('A','A','A','A','A','A','A','A','A','A','a','a','a','a','a','a','a','a','a','a') out=c()
while(sum(gen=='A') > 0 && sum(gen=='a') > 0) { Så længe der både er A-er og a-er...
out = c(out,sum(gen=='A')) ... tæl A-er og husk antallet ...
gen = sample(gen,replace= TRUE) ... og lav en ny generation.
}
R - indledende eksempler
plot(out,type='l',col='blue',ylim=c(0,20))
0 5 10 15 20 25 30
05101520
Index
out
I fraværet af mutationer, vil A = 0 eller a = 0 være
absorberende tilstande.
R - indledende eksempler
Ved at tænke på et valg af A som “succes” ses at antallet af A-er binomialfordelt, hvor sandsynlighedsparameteren er frekvensen af A-er i den foregående generation.
N = 2000 frekvens = 0.5 out=c(frekvens)
while (frekvens > 0 && frekvens < 1) {
antalAer = rbinom(1,N,frekvens) rbinom simulerer binomialfordelingen.
frekvens = antalAer/N beregn ny frekvens out = c(out,frekvens)
}
R - indledende eksempler
plot(out,type="l",col='blue',ylim=c(0,1))
0 500 1000 1500 2000 2500
0.00.20.40.60.81.0
out
Estimation af absorptionssandsynligheden
Med et simpelt R program kan vi estimere sandsynligheden for absorption i 1:
N = 2000 R = 10000
out = rep(0,R) Initialiserer en vektor med R nuller
for (i in 1:R) { frekvens = 0.5
while (frekvens > 0 && frekvens < 1) { frekvens = rbinom(1,N,frekvens)/N }
out[i] = frekvens }
mean(out) mean beregner gennemsnittet
Output:
[1] 0.5033
Estimation af ventetiden på absorption
Et lidt mere kompliceret spørgsmål er følgende: Hvordan sammenhængen mellem middel-ventetiden på absorption og frekvensen?
Spørgsmålet kan belyses vha. nedenstående stykke kode:
for (initial_frekvens in 1:9) { Gennemløb 1 . . .9 absorptions_tider = rep(0,R)
for (i in 1:R) {
c = 0 tæller til antal gennemkørsler af løkken
frekvens = initial_frekvens/10 start-frekvensen while (frekvens > 0 && frekvens < 1) {
frekvens = rbinom(1,N,frekvens)/N c = c + 1
}
absorptions_tider[i] = c }
out[initial_frekvens] = mean(absorptions_tider)
}
Estimation af ventetiden på absorption
Plot af resultatet:
0.0 0.2 0.4 0.6 0.8 1.0
5001000150020002500
x
theoabstime(x)
●
●
●
●
●
●
●
●
●
Den indtegnede kurve er grafen for funktionen
p 7→ −2N(p log(p) + (1 − p) log(1 − p)) .
Mutationer indføres
Vi antager nu at der i hver generation er sandsynlighed p A→a
for at mutere fra A til a, og p a→A for at mutere fra a til A. Antallet af A-er findes nu som
#A-er = #A-er valgt i foregående generation, som ikke muterede
+ #a-er valgt i foregående generation, som muterede. +(N − #A-er) valgt i foregående generation, som muterede.
Idet f betegner frekvensen af A-er i foregående generation, resulterer dette i en binomialfordeling med
sandsynlighedsparameter:
f(1 − p A→a ) + (1 − f )p a→A .
Når mutationer kan forekomme vil A = 0 og a = 0 ikke længere
være absorberende tilstande.
Histogrammer med mutationer
N = 200; T = 300; repetitioner = 50000 T=300generationer p_A_to_a = p_a_to_A = 0.002 Dette erpA→a ogpa→A
out =matrix(0,repetitioner,T) out[,1] = N/2
for(i in 1:repetitioner) { for (j in 2:T) {
out[i,j] =rbinom(1,N,(out[i,j-1]/N)*(1-p_A_to_a)+(1-out[i,j-1]/N)*p_a_to_A) }
}
hist(out[,10]/N) Histogram for den 10ende søjle, dvs. den10ende generation.
Generation = 10
Frekvens af A−er
Tællinger
0.0 0.2 0.4 0.6 0.8 1.0
05001000150020002500
Generation = 50
Frekvens af A−er
Tællinger
0.0 0.2 0.4 0.6 0.8 1.0
020040060080010001200
Generation = 300
Frekvens af A−er
Tællinger
0.0 0.2 0.4 0.6 0.8 1.0
0200400600
Histogrammer med mutationer
Histogrammer som funktion af tiden.
Histogrammer med mutationer
Histogrammer, alternativ startbetingelse.
To populationer
Vi forestiller os nu a vores oprindelige population af størrelse N splitter op i to populationer af størrelse N /2 hver.
Vi indfører derfor migrationssandsynligheder p 1→2 og p 2→1 . For en population fandt vi en binomialfordeling med
sandsynlighedsparameter
f(1 − p A→a ) + (1 − f )p a→A .
Med to populationer bliver antal A-er i population 1 igen binomialfordelt, nu med sandsynlighedsparameter
f 1 (1 − p A→a )(1 − p 1→2 ) + (1 − f 1 )p a→A (1 − p 1→2 )+
f 2 (1 − p A→a )p 2→1 + (1 − f 2 )p a→A p 2→1 og vi får et analogt udtryk for den tilsvarende
sandsynlighedsparameter for population 2.
Histogrammer med mutationer - to populationer
Histogrammer for to populationer.
Heterozygositet
Vi indfører nu heterozygositeten til at være sandsynligheden for at udtrække 2 forskellige alleler fra populationen, dvs.
H n = 2(# A-er)(N − # A-er)
N 2 ,
og ser på hvordan denne opfører sig som funktion af tiden.
Histogrammer og estimeret heterozygositet
Efterskrift
I Wright-Fisher modellen er udgangspunktet nutiden, og modellen beskriver den fremtidige udvikling.
I forbindelse med inferens vil vi i praksis betragte en model,
hvor tiden “kører” baglæns, dvs. hvor modellen beskriver
udviklingen bagud i tid.
“Baglæns tid”
1 2 3 4 5 6 7
Fortid
Nutid
Afslutning
Scally et al. Nature, 2012.
Andersen et al. Journal of Mathematical Biology (2014)