Skattning av avstånd mellan arter i fylogenetiska träd

U.U.D.M. Project Report 2018:28

Examensarbete i matematik, 15 hp Handledare: Ingemar Kaj

Examinator: Martin Herschend Juni 2018

Department of Mathematics

Skattning av avstånd mellan arter i fylogenetiska träd

Linnéa Eriksson

Sammanfattning

I det här arbetet beräknas och visas det om skattning av avståndet mellan gen-

frekvenser. De matematiska modellerna kan tillämpas inom fylogeni. Modellerna

som arbetet tar upp är JC69, Jukes-Cantor, och K80-modellen, Kimura. De två

modellerna studeras steg för steg och tillämpas därefter på en människas genfre-

kvens mot fem olika djur. De fem djuren är schimpans, gorilla, bonobo, gibbon

och ett utstickande djuret som är lejon. Genfrekvenserna från djuren som är

jämförda med människans är hämtade ifrån GeneBank. Genfrekvenserna stude-

ras och data tillämpas sedan på de matematiska modellerna. Beräkningar och

grafer har utförts i datorprogrammet MatLab. Slutligen så jämförs alla beräk-

ningar med varandra och de diskuteras hur man skulle kunna gå tillväga för att

utveckla arbetet.

Innehåll

1 Inledning 4

2 Biologisk bakgrund 5

2.1 Fylogeni . . . . 5

2.2 DNA . . . . 5

2.3 Protein och nukleotider . . . . 6

2.4 Kodonbias . . . . 7

3 Modeller för nukleotidsubstitution 7 3.1 JC69 (Jukes and Cantor 1969) . . . . 8

3.2 K80 (Kimura 1980) . . . 11

3.3 Generellt för båda modellerna . . . 14

3.4 Avståndsuppskattning med UNREST . . . 15

4 Maximum likelihood-metoden 15 4.1 JC69 . . . 16

4.2 K80 . . . 17

5 Uppbyggnad av fylogenetiska träd 17 5.1 Avstånd mellan arter . . . 18

5.1.1 Minstakvadratmetoden . . . 18

5.2 Maximum likelihood-metoden - er generationer . . . 19

5.2.1 Likelihood beräkningar på träd . . . 19

6 Resultat 20 6.1 Människa (Homo sapiens) D38112 mot Schimpans (Pan troglo- dytes) . . . 20

6.1.1 JC69 modellen . . . 21

6.1.2 K80-modellen . . . 23

6.2 Människa (Homo sapiens) D38112 mot Gorilla (Gorilla gorilla) . 25 6.2.1 JC69 modellen . . . 25

6.2.2 K80-modellen . . . 27

6.3 Människa (Homo sapiens) D38112 mot Bonobo (Pan paniscus) . 29 6.3.1 JC69 modellen . . . 30

6.3.2 K80-modellen . . . 32

6.4 Människa (Homo sapiens) D38112 mot Svarthandad Gibbon (Hy- lobates agilis) . . . 34

6.4.1 JC69 modellen . . . 34

6.4.2 K80-modellen . . . 36

6.5 Människa (Homo sapiens) D38112 mot Lejon (Panthera leo) . . . 38

6.5.1 JC69 modellen . . . 39

6.5.2 K80-modellen . . . 41

6.6 Jämförelse av resultat . . . 43

7 Diskussion 45

8 Bilagor 46 8.1 Tabell för 95%-kondensintervall - normalfördelnings kvantiler . 46

8.2 Tabell för χ

. . . 46

8.3 Genfrekvenser från GeneBank . . . 47

8.3.1 Människa (Homo sapiens) D38112 mot Schimpans (Pan troglodytes troglodytes) . . . 47

8.3.2 Människa (Homo sapiens) D38112 mot Gorilla (Gorilla gorilla) . . . 48

8.3.3 Människa (Homo sapiens) D38112 mot Bonobo (Pan pa- niscus) . . . 49

8.3.4 Människa (Homo sapiens) D38112 mot Svarthandad Gib- bon (Hylobates agilis) . . . 50

8.3.5 Människa (Homo sapiens) D38112 mot Lejon (Panthera leo) 51 8.4 Matlab-kod för uträkningar för modellerna . . . 52

8.4.1 Matris - Människa mot Schimpans . . . 52

8.4.2 Matris - Människa mot Gorilla . . . 52

8.4.3 Matris - Människa mot Bonobo . . . 53

8.4.4 Matris - Människa mot Svarthandad Gibbon . . . 53

8.4.5 Matris - Människa mot Lejon . . . 54

8.4.6 Kod för JC69 . . . 54

8.4.7 Kod för K80 . . . 56

9 Referenser 58 9.1 Referenser för matematiska modeller . . . 58

9.2 Referenser för fakta . . . 58

9.3 Referenser genfrekvens . . . 58

9.4 Referenser för bilder . . . 58

1 Inledning

Att beräkna avståndet mellan två genfrekvenser är en relativt enkel fylogenetisk analys men ändå väldigt viktig. En viktig del är beräkningar av avstånd mellan sekvenspar. Vilket är de första steget i metoden för konstruktion av avstånds- matrisen inom fylogeni. De består av att klusteralgoritmer som konverterar en avståndsmatris till ett fylogenetiskt träd. De andra viktiga är modeller för mar- kovprocesser av nukleotidsubstitution. Det används i avståndsberäkningarna från basen av maximum likelihood och bayesiansk analys av multipla sekvenser i fylogeni.

I det här arbetet har jag till stor del utgått från boken Computational Mo- lecular Evolution av Ziheng Yang. Det är en modern bok som bygger på statis- tiska och beräkningsmässiga metoder som används i molekylär evolutionsanalys, såsom maximum likelihood, markovprocesser och bayesianska statistik. I boken analyseras molekylär sekvensdata och som vi under de senaste åren fått extremt mycket mer förståelse för. Boken går inte in på djupare matematiska bevis utan håller sig till metoder och hur de beräknas. Det nämns även lite om hur man går till väga för att påbörja byggandet av ett fylogenetiskt träd.

Metoderna som studeras i det här arbetet är de två modellerna JC69 och

K80. Det används olika metoder så som markovprocesser, maximum likelihood-

metoden och avståndsmetoden. De olika genfrekvenserna som studerats i ar-

betet kommer från GeneBank och beräkningarna som genomförts har gjorts i

datorprogrammet MatLab.

2 Biologisk bakgrund

Genom matematisk statistik analyseras biologisk data för att få fram sannolik- heten av önskad hypotes. För att sedan studera den samt se vad resultatet visar och om hypotesen kan förkastas eller inte.

2.1 Fylogeni

Fylogeni är en studie om organismers släktskap där resultaten sammanställs med fylogenetiska träd. Idag studeras släktskap mellan organismer genom att jämfö- ra deras DNA. Längre tillbaka när människan inte hade någon större vetskap om DNA, studerade man de olika organismernas yttre och morfologiska egenskaper.

Under de senaste åren har vetenskapen om DNA och den molekylära evolutio- nen ökat explosionsartat. Detta då kunskapen har ökat något enormt inom det tekniska, vilket gör att det nu går mycket snabbare att ackumulera genetisk sekvensdata, vår förbättring inom hårdvara och mjukvara samt utvecklingen av analysmetoder.

Den stora ökningen av genomisk data kräver kraftfulla statistiska modeller och datorer för att de ska kunna analyseras och tolkas. Tre termer som ofta används inom fylogeni är monofyli, parafyli och polyfyli. Monofyli är de som omfattar ättlingar, det vill säga de närmsta individerna med gemensam stamfa- der och gemensamma förfäder. Parafyli är när en grupp bestående av ättlingar till en stamfader men i denna grupp ingår inte alla ättlingar utan vissa kan ute- slutas på grund utav olika anledningar. Polyfyli är en grupp som är besläktade men inte nära, det vill säga de har en avlägsen gemensam stamfader.

De vetenskapliga metoderna som används inom fylogenetik brukar gruppe- ras i vad som benämns kladistik. Skillnaden mellan fylogenetik och kladistik är att fylogenetik kan innehålla hypoteser om släktskap, medan kladistik istället tillämpas mer vetenskapligt som till exempel i matematiska modeller. Vanli- ga vetenskapliga modeller som används inom kladistik är maximum likelihood- metoden och markovprocesser med en bayesiansk inferens. Markovprocesser och maximum likelihood-metoden är de metoder som kommer att studeras i det här arbetet.

Inom fylogeni och kladistik studeras homologa egenskaper hos organismer.

Organismer som har homologa egenskaper anses vara närmare besläktade och tvärtom, färre likheter mer avlägsna från varandra. Inom detta konstruerar man träd för att på ett enkelt sätt se hur organismer är besläktade, de benämns fylogenetiska träd eller kladogram. Dessa två är väldigt lika varandra, de som bland annat skiljer dem åt är att i ett fylogenetiskt träd indikerar grenarna på olika tidsförhållanden.

2.2 DNA

DNA är en förkortning av deoxyribonucleic acid och är det ämne i en organism

som bär på den genetiska informationen. DNA-molekylens viktigaste funktion är

att lagra information om organismens funktioner och utveckling. DNA innehåller

all information om hur organismen ska konstruera och hur den ska bygga upp alla ämnen. Därför kan DNA kallas för kroppens alldeles egna receptbok. En DNA-molekyl har två strängar, så kallade polymer, som i sin tur är uppbyggda av nukleotider. En nukleotid består av en kvävebas och en pentos. Där det är kvävebaserna som innehåller den genetiska koden. Det nns fyra olika typer av kvävebaser, Adenin (A), Cytosin (C), Guanin (G) och Tymin (T). De fyra nukleotider kan inte kopplas samma hur som helst, adenin och tymin kopplas ihop samt cytosin och guanin kopplas ihop.

Fig.1 En DNA spiral som delas och som visar alla fyra nukleotider samt hur de kopplas samman.

2.3 Protein och nukleotider

Ett protein är en lång kedja bestående av aminosyror. En aminosyra är kemis- ka föreningar mellan en aminogrupp och en karboxylgrupp. Det nns en stor mängd olika aminosyror men alla nns inte levande i organismers celler. När det pratas om levande organismer säger man att det existerar 20 aminosyror, det

nns dock enstaka undantag för några få organismer. Ett protein byggs upp in- uti en cell i två steg. I det första steget transkriberas proteinet och det benämns även för RNA-syntes. Detta är en process där den genetiska informationen i en cells DNA kopieras och skapar ett RNA. RNA är som DNA uppbyggt av nuk- leotider. Nukleotiderna är nästintill lika som de som används för DNA, det som skiljer är att tymin (T) har ersatts av uracil (U). När själva transkriptionen sker delar sig DNA-strängen för att den ska kunna bilda en mall för RNA:t.

Nukleotiderna A, C, G och T i DNA-kedjan kommer att ge upphov till U, G, C och A på motsvarande plats i RNA-molekylen. Denna typ av RNA benämns för mRNA som är förkortning för messenger RNA och agerar som förnamnet säger som budbärare mellan cellkärnan och ribosomerna. I andra steget translaterar mRNA i ribosomerna till aminosyror. Det är alltså här som aminosyrorna sätts samman till det färdiga proteinet. Translationen sker så att nukleotiderna kopp- las ihop och läses av tre och tre, där tre nukleotider kodar tillsammans för en specik aminosyra. En grupp av tre nukleotider benämns för ett kodon och det

nns 4

= 64 möjliga kombinationer. Det nns dock bara 20 olika aminosyror,

så olika kodon kan koda för samma aminosyra. Det nns ett startkodon som

startar translationen och tre stoppkodon som gör att translationen avbryts.

2.4 Kodonbias

Det nns fördelar med att en aminosyra svara för er än ett kodon, en anledning är att den blir mer tålig mot mutationer. Det förekommer nämligen främst att det är den sista nukleotiden som har översatts felaktigt utav de tre nukleotider- na. Organismer fungerar så nurligt att den sista nukleotiden sällan har någon större betydelse för vilken aminosyra den kodar för. Det är alltså de två första nukleotiderna som till största del är avgörande för vilken aminosyra den kodar för. Studeras till exempel aminosyran alanin och dess genetiska kod så är den GCU, GCC, GCA och GCG. Här ser man tydligt att de två första nukleotider- na är densamma varav den sista varierar och alla kodar ändå för alanin. Det är dock inte alltid så simpelt, organismer har listigt nog gjort att aminosyror med liknande kodon generellt har relativt lika egenskaper. Vilket därför sällan gör någon större skillnad om ett kodon blir fel då aminosyran oftast har likvärdiga funktioner som den tilltänkta aminosyran. Därav kan proteinet relativt ofta fun- gera som de ska ändå. Kodonbias är alltså de praktiska som organismen skapat som gör att relativt små skillnader i ett kodon sällan har någon större betydelse för de kodande DNA:t. Kodonbias kan även förekomma i övergångarna mellan olika nukleotider och göra att det inte sker likformigt. Övergångarna delar man i transition och transversion. Först delas nukleotiderna in i två grupper, puriner och pyrimidin. Nukleotiderna A och G är puriner som är heterocykliska kväve- föreningar som är uppbyggda av två ringar. C och T är pyrimidin och är även dem heterocykliska föreningar men är endast uppbyggd av en ring istället. Över- gångarna inom grupperna purin och pyrimidin är transitioner och övergångar mellan de två grupperna är transversioner.

Fig.2: Bilden visar vilka övergångar mellan nukleotider samt vilka som är transition och transversion. De blåa pilarna är transversion och de röda är transition.

3 Modeller för nukleotidsubstitution

Här kommer två modeller för nukleotidsubstitution att studeras steg för steg.

Den första heter Jukes-Cantor modellen, JC69, och är en av enklare modeller för

nukleotidsubstitution. Den andra heter Kimura modellen, K80, och är fortfaran- de relativt enkel men något mer avancerad än Jukes-Cantor modellen. Figur 3, nedan, illustrerar övergångarna mellan de fyra nukleotiderna samt hur transition och transversion har lite olika betydelse mellan de två modellerna. För JC69- modellen har samtliga övergångar samma frekvens, a. För K80-modellen skiljer sig övergångarna, transitioner har frekvens a och transversioner har frekvens b.

Fig.3: Visar hur övergångarna är för modell K80. JC69 är liknande bara att alla övergångar är α, det vill säga β = α.

3.1 JC69 (Jukes and Cantor 1969)

JC69 antar att alla nukleotidsekvenser har samma frekvens, λ, av en förändring till en annan nukleotid. Frekvensen q

= ¨ ogonblicksf rekvensen av substitution f r˚ an nukleotid i till j, d¨ ar i, j = T, C, A och G. Matris (1), nedan, har ordningen T, C, A och G för nukleotiderna. Varje matrisrad måste ha summan noll. Den totala substitutionskvoten för bytet av nukleotid, i, är 3λ vilket i matrisen står för q

. Det är −q

som motsvarar substitutionfrekvensen för nuk- leotid, i, det vill säga frekvens det tar för markovkedjan att lämna tillståndet i.

Frekvensmatrisen är

Q = {q

} =









−3λ λ λ λ

λ −3λ λ λ

λ λ −3λ λ

λ λ λ −3λ









(1)

Övergångsmatrisen är P (t) = {p

(t)} . Övergångssannolikheten, p

(t) , är san- nolikheten där given nukleotid, i, vill bli nukleotid, j, över tiden, t. Beräkningen på övergångsmatrisen ger som följer

P (t) = { p

(t) } = e

=









p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t)







 (2)

där övergångssannolikheten är

( p

(t) =

+

e

p

(t) =

−

e

(3) Beräkningar görs på övergångsmatrisen P (t) och en matris exponentiellt genom taylorutveckling. Denitionen för taylorutveckling följer

P (t) = e

= I + Qt + 1

2! (Qt)

+ 1

3! (Qt)

+ 1

4! (Qt)

+ . . . (4) En taylorutveckling på nukleotider är varken avancerad eller tidskrävande då matrisen generellt är förhållandevis liten. Emellertid kan denna metod bli mer kostsam och ostabil om man gör det för någon aminosyra eller för ett kodon då det ger en mycket större matris. En matris för en aminosyra är storlek 20 x 20 och för ett kodon 61 x 61. Från matris (2) kommer i att för varje plats vara någon nukleotid i en lång sekvens under tiden, t. Den andra nukleotiden j i en sekvens kommer att vara p

(t) , där j = T, C, A, G. Summan av varje radmatrisen är ett, P (t) = 1. För tiden noll, t = 0, är blir övergångsmatrisen identitetsma- trisen, P (0) = I. De nns två generella modeller för markovkedjor. Den första är den generella tidsövergångsmodellen och den andra är den generella otvung- na modellen. Notera att markovprocesser klassiceras beroende på om tiden på tillståndet är diskret eller kontinuerligt. Den teori som visas här, för JC69, är en relativt enkel modell där utbytet sker mellan aminosyror och kodon. När t → ∞ är p

(t) =

för alla i och j. Detta visar när en substitution har inträat många gånger på varje plats, så att den slutliga nukleotiden är slumpmässig med sannolikheten

för varje nukleotid oberoende från vart man började. Sannolik- heten att kedjan är i tillståndet j när t → ∞ betecknas för π

. Distributionen är (π

, π

, π

, π

) och benämns för limiting distribution. För JC69 är π

=

för varje nukleotid j, där jämviktdistributionen blir π = (

,

,

,

) . Detta ger πQ = 0 givet att P

π

= 1 . Om det nns en markovkedja med era tillstånd används följande ekvation, även kallad för Chapman-Kolmogorov teoremet

p

(t

+ t

) = X

k

p

(t

)p

(t

) (5) Sannolikheten att nukleotid i blir nukleotid j under tiden t

+ t

är summan av alla möjliga tillstånd, k, vid varje mellanliggande tidpunkt t

. Det är avståndet mellan dessa två sekvenser som ska beräknas. Från frekvensmatrisen, matris (1), får man den totala substitutionsfrekvensen för någon nukleotid, som är 3λ . Därav kan avståndet mellan två sekvenser beräknas till d = 3λt. Där d är avståndet, t är tiden och λ är frekvensen. Antag att x utav n platser är olika mellan två sekvenser, då kommer proportionen av dierensen av platserna att bli ˆp =

. Detta är sannolikheten, p, för att en plats har olika nukleotider mellan de två sekvenserna med ett avstånd, d, som ger följande

p = 3p

(t) = 3 4 − 3

4 e

=  d 3 = λt



= 3 4 − 3

4 e

(6)

Beräkning för att räkna ut den uppskattade avståndet är

ˆ p − 3

4 = 3

4 e

⇔ 1 − 4

3 p = e ˆ

Vidare förenklingar ger

log(1 − 4

3 p) = log(e ˆ

) ⇔ log(1 − 4

3 p) = − ˆ 4 ˆ d 3

− 3

4 log(1 − 4 3 p) = ˆ ˆ d

Följande blir den slutgiltiga uppskattningen för avståndet d = − ˆ 3

4 log(1 − 4

3 p) ˆ (7)

När ˆp >

går de skattade avståndet inte att tillämpas, två slumpmässiga sekvenser bör alltså ha omkring 75% olika platser. När ˆp <

är de skattade avståndet oändligt. Sannolikheten, p, är binomial i förhållande till variansen,

ˆ p(1− ˆp)

n

. Variansen av de skattade avståndet, ˆ d , ska nu härledas, där ˆ d är en funktion av den skattade sannolikheten, ˆp. Gauss-approximationen används för att räkna ut variansen.

var( ˆ d) = var(ˆ p) · | d ˆ

dˆ p |= p(1 − ˆ ˆ p)

n · 1

(1 −

)

(8) Gauss-approximationen används som en generell riktlinje för att derivera vänte- värdet, variansen och kovariansen av funktion med slumpmässiga variabler. Där en icke-linjär funktion, f(x), där x är en slumpmässig variabel som har vänte- värdet är µ och variansen är σ

. Följande gäller E(f(x)) 6= fE((x)). När n är ett positivt heltal, kan taylorutvecklingen skrivas som följande med ordningen n och funktionen, f,

f (x) = T

(x) + R

(x) där

T

(z) = f (a) + f

(a)

1! (z − a) + f

(a)

2! (z − a)

+ . . . + f

(a)

n! (z − a)

(9) Taylorutvecklingen ska nu tillämpas på Gauss-approximation. Taylorutveckling av f(x) runt väntevärdet µ ger

f = f (x) = f (µ) + df (µ)

dx (x − µ) + d

f (µ)

2! dx

(x − µ)

+ . . . (10)

Funktionen, f, och derivatorna är ekvivalent med x = µ. Alla termer med

exponent tre eller högre ger ett väntevärde för funktionen. Det approximerade

väntevärdet för funktionen, f, blir

E(f ) ≈ f (µ) + 1 2

d

f (µ) dx

σ

Där E(x − µ) = 0 och E(x − µ)

= σ

. Derivatan är ekvivalent med x = µ och de är konstant där när man tar förväntade värden över x. Den approximerande variansen av funktionen, f, och den uppskattade parametern x.

var(f ) ≈ E(f − E(f ))

≈ σ

 df (µ) dx



Efter att ha räknat ut variansen av den skattade sannolikheten, ˆp, var(ˆp) =

ˆ p(1− ˆp)

n

, och variansen av det skattade avståndet, ˆ d , var( ˆ d) =

·

, kan man beräkna derivatan av dem,

=

(1−^{4 ˆ}₃^p)

. Slutligen tillämpas detta tillsammans med ett approximerat 95%-kondensintervall, ˆ d ± λ

. ε . Där ε är de standard felet, ε = q

var( ˆ d) , och där signikansnivån på 95% ger λ

= 1.96 . Olika signikanta nivåer och dess värden kommer från tabell som ligger under bilagor, bilaga 8.1.

3.2 K80 (Kimura 1980)

I K80 modellen nns en substitution mellan antingen två pyrimidin (hetero- cyklisk förening, cytosin och tymin) T ↔ C eller mellan två puriner (hetero- cyklisk kväveförening, har två ringar, adenin och guanin) A ↔ G. När någon av dessa två sker benämns det för transition. Substitutioner sker mellan pyrimidin och puriner (T, C ↔ A, G) och det benämns för transversioner. I verkligheten uppkommer transitioner med högre frekvens än transversioner. Alltså notera att transitionen och transversionen inom biologin inte har exakt samma sannolikhet som för modellerna. Substitutionsfrekvensen för transitionen kallas för α och för transversionen β. Frekvensmatrisen blir som följer

Q = {q

} =









−(α + 2β) α β β

α −(α + 2β) α β

β α −(α + 2β) α

β β α −(α + 2β)







 (11) Den totala substitutionsfrekvensen för någon nukleotid är α+2β. Där avståndet mellan två sekvenser multipliceras med tiden t, vilket ger avståndet d = (α + 2β)t . Där αt är det förväntade värdet transitioner per plats och 2βt är det förväntade värdet för transversioner per plast. Oftast används avståndet, d, eller transitions-/transversionsfrekvenskvoten, κ =

. Jämviktsfördelningen för K80 är identisk som för JC69, alltså π = (

,

,

,

) . Där πQ = 0 givet att P

i

π

= 1 . Övergångsmatrisen är följande

P (t) = {p

(t)} = e

=









p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t)







 (12)

De tre olika övergångssannolikheterna, p, som nns i matrisen beräknas via taylorutveckling, se ekvation (9), och blir med de nya värdena.



 

 

p

(t) =

+

e

+

e

=

+

e

+

e

p

(t) =

+

e

−

e

=

+

e

−

e

p

(t) =

−

e

=

−

e

Summan av en radmatris måste bli värdet ett, det vill säga p

(t) + p

(t) + (13) 2p

(t) = 1 . Denna datasekvens kan nu delas in i andelar av transitional och transversional dierens, de kommer att få betäckningarna S och V. Genom symmetrin i modellen och matris (12) blir sannolikheten för uppkomsten av en plats nukleotider den transitionala dierensen E(S) = p

(t) och transversionala dierensen E(V ) = 2p

(t). Där det skattade avståndet, ˆ d , samt de skattade transitions-/transversionsfrekvenskvoten, ˆκ, ger

p = p

(t) + 2p

(t) = 1 4 + 1

4 e

− 1

2 e

+ 2( 1 4 − 1

4 e

)

= 1 4 + 1

4 e

− 1

2 e

+2( 1 4 − 1

4 e

) = 3 4 − 1

4 e

− 1

2 e

Det skattade avståndet, ˆ d , blir d = − ˆ 1

2 log(1 − 2S − V ) − 1

4 log(1 − 2V ) (14)

Där det skattade transitions-/transversionsfrekvenskvoten, ˆκ, blir

ˆ

κ = 2log(1 − 2S − V )

log(1 − 2V ) − 1 (15)

Transitionsavståndet är ekvivalent med αt och transversionsavståndet med 2βt och är skattade till följande

αt = − b 1

2 log(1 − 2S − V ) + 1

4 log(1 − 2V ) (16)

2 b βt = − 1

2 log(1 − 2V ) (17)

Transitionsavståndet gäller endast om 1−2S −V > 0 samt 1−2V > 0. S och V får följande varianser var(S) =

och var(V ) =

. Det ger en kovari- ans på cov(S, V ) = −

. Därefter används Gauss-approximationen, se ekvation (10), och deriverar varians-kovariansmatrisen. Varians-kovariansmatrisen ser ut som följer

var( S V ) =

S(1−S) n

−SV n

−

!

(18)

där n står för antal platser i sekvensen. En skattning görs av ˆ d och ˆκ på en funktion av S och V . Tillämpningen blir följande

var( d ˆ ˆ

κ ) = J · var( S

V ) · J

(19)

J står för en Jacobimatris av en anpassad storlek m x n. Här blir Jacobianen följande

J =

∂ ˆd

∂S

∂ ˆd

∂ ˆκ ∂V

∂S

∂ ˆκ

∂V

!

=

1 1−2S−V

1

2(1−2V )

+

−

−

+

!

(20) Det gör att att variansen av ˆ d slutligen kan deriveras. Så

var(f ) ≈

n

X

i=1 n

X

j=1

cov(x

, x

)( ∂f

∂x

)( ∂f

∂x

) (21)

som är variansen av ett enkelvärdesfunktion av f(x) approximerat av x. Där cov(x

, x

) är kovariansen av x

och när i 6= j och när i = j blir det variansen istället. Därefter får man

var( ˆ d) =( ∂ ˆ d

∂S )var(S) + 2 · ∂ ˆ d

∂S · ∂ ˆ d

∂V · cov(S, V ) + ( ∂ ˆ d

∂V )

var(V )

= a

S + b

V − (aS + bV )

 /n

(22)

där a och b står för följande

a = (1 − 2S − V )

(23)

b = 1

2 (1 − 2S − V )

+ (1 − 2V )

 (24)

Till sist kan man använda det man beräknat och tillämpa de i ett approximerat 95%-kondensintervall ˆ d ± λ

. ε . Avståndet och standard felet räknas ut ifrån S och V . Det går även att studera transitionsfrekvensen för att se hur mycket större eller mindre den är jämfört med transversionsfrekvensen.

3.3 Generellt för båda modellerna

Låt tillståndet i kedjan vid tidpunkt t så den blir X(t). X(t) är en av de fy- ra nukleotiderna A,T,C eller G. Antag att alla positioner i en DNA-sekvens utvecklas oberoende och att markovprocessen används för att beskriva nukleo- tidsubstitutionerna för alla positioner. Där är P r {X(t + t | X(t) = i} = q

t om markovprocessen q

beror av tiden t. Beror q

inte av tiden, t, säger man att den är tidshomogen. Den generella modellen utan några begränsningar av uppbyggnaden utav frekvensmatrisen, Q, kommer att bestå av 12 fria paramet- rar. Den angivna frekvensmatrisen, Q, över någon tid t > 0 : P (t) = {p

(t)} , där p

(t) = P r {X(t) = j | X(0) = i} . P (t) ger alltså följande ekvation

dP (t)

dt = P (t)Q, (25)

med ett randvillkor på P (0) = I, där I är identitetsmatrisen. Detta ger i sin tur lösningen

P (t) = e

(26)

Frekvensmatrisen, Q, och tiden, t, är en produkt där Q varierar i olika skal- faktorer att den genomsnittliga frekvensen blir ett. Markovkedjan X(t) har en initial fördelning π

= (π

, π

, π

, π

) , medan tiden t har fördelningen π

= (π

, π

, π

, π

) , vilket ger följande

π

= π

P (t) (27)

Ett exempel om man tar ekvation (27) med nukleotiden T som slutpunkt och ett initialvärde på noll. Får man ekvationen π

= (π

, π

, π

, π

) . När den initiala och slutliga fördelningen är ekvivalenta, det vill säga π

= π

, så kommer kedjan att stanna i fördelningen i en oändlighet. Då säger man att kedjan är stationär eller att den är i jämvikt. Man säger även att fördelningen π är stationär eller i steady-state fördelning. Markovkedjan gör att alla tillstånd kan anta vilket annat tillstånd inom en ändlig tid med en positiv sannolikhet.

Denna kedja säger man är irreducibel och har då en unik stationär fördelning, vilket också är begränsad fördelning när tiden t → ∞. Från ekvation (27) är följande ekvivalent

πQ = 0 (28)

Notera att de totala ödet av något j är P

π

q

medan det totala utödet är för något j är −π

q

. När kedjan är stationär kommer detta tillsammans med ekvation (28) att vara identiska, alltså P

π

q

= 0 för något j. Ekvation (28) tillsammans med π

> 0 och P

j

π

= 1 ger oss möjlighet att bestämma

den stationära fördelningen från Q för någon markovkedja.

3.4 Avståndsuppskattning med UNREST

UNREST är en generell modell för nukleotidsubstitution med en frekvensmatris, Q , utan några större begränsningar och med 12 parametrar. Frekvensmatrisen, Q , denieras av den relativa frekvensen där 11 parametrar är involverade. Model- len implementerades av Yang (1994b) för att uppskatta avståndet av sekvenser som använder två grenlängder, t

och t

. Maximum likelihood-metoden ger den multinomiella sannolikheten med 16 olika cellerna, där de 16 cellerna motsvarar 16 möjliga kombinationer. Låt funktionen f

(t

, t

) vara sannolikheten för den ij -te cellen, det vill säga den sannolikhet för att någon plats har nukleotid i i ena sekvens och nukleotid j i den andra sekvens. Då de fyra möjliga nukleotiderna härstammar från förfäder måste den genomsnittliga beräknas över dem

f

(t

, t

) = X

k

π

p

(t

)p

(t

) (29) Låt n

vara antalet platser i den ij-te cellen. Då blir logaritmen av maximum likelihood-metoden som följer

`(t

, t

, Q) = X

i,j

n

log{f

(t

, t

)} (30) Frekvensparametrarna π

, π

, π

, π

denieras från frekvensmatrisen, Q, med hjälp av ekvation (30) och de är inte fria parametrar. Det nns dock två problem med denna modell som därför inte alltid gör att den är helt lämplig att tillämpa. Det ena problemet är att den numeriska metoden som används för att hitta maximum likelihood-metodens parametrarna där ingen analytisk lösning verkar möjlig. Egenvärdena för frekvensmatrisen Q tar nämligen inte hänsyn till komplexa tal. Den andra anledningen är att den typiska datamängderna sällan är tillräckliga för att ge otillräckligt med information för att kunna skatta parametrarna.

4 Maximum likelihood-metoden

Generellt används maximum likelihood som en metod för att skatta paramet- rar i en modell och för att testa hypoteser om parametern. Denna metod har många användningsområden och inom molekylärfylogeni har den en viktig roll.

Här används maximum likelihood-metoden för att uppsatta avståndet i en se-

kvens. Låt X vara vår data och θ den parameter som man vill skatta. Maximum

likelihood-funktionen betecknas som följande L(θ; X) = f(θ | X). Den kan med

ord förklaras så att sannolikheten av den observerade informationen är X med

en studerad funktion av en okänd parameter θ, med en given data. Likelihood-

principen säger att maximum likelihood-funktionen har all information i data

om θ.

4.1 JC69

Vid användning av maximum likelihood-metoden för JC69 skattas avståndet mellan sekvenser och parametrar, där avståndet är d. Data för två sekvenser som vardera har n platser och x antal skillnader mellan de två sekvenserna.

Detta är sannolikheten p för att en plats har olika nukleotider mellan de två sekvenserna med en avståndet d, som ger följande

p = 3p

(t) = 3 4 − 3

4 e

(31)

Sannolikheten för den observerade data, x, som är antal skillnader mellan de två sekvenserna och består av n antal platser, får man genom den binomial sannolikheten.

L(d; x) = f (x | d) = Cp

(1 − p)

= C( 3 4 − 3

4 e

)

( 1 4 + 3

4 e

)

(32) Sannolikheten av den observerade data, x, skrivs som en funktion utav parame- tern d, det vill säga avståndet. Värden för avståndet, d, med ett högt värde av maximum likelihood, L, stöds bättre än för låga värden för maximum likelihood, L . Ekvation (32) ska nu kompletteras. Först adderas den binomiala koecien- ten, C = h

n!

x!(n−x)!

i, men eftersom det är en konstant och kan den förkastas.

Samma denition kommer därför att användas för samtliga substitutionsmodel- ler och där det nns 16 möjligheter istället för två stycken som i ekvation (32), p och 1 − p. I JC69 nns fyra konstanta mönster (TT, CC, AA, GG) där alla har samma sannolikhet att inträa, där det är lika för de andra 12 möjligheterna (TC, TA, TG etc.). Detta är en omdenierad multinomial sannolikhet för 16 celler

L(d; x) = ( 1 4 p

)

( 1

4 p

)

= ( 1 16 − 1

16 e

)

( 1 16 − 3

16 e

)

(33) Parametrarna p

och p

kommer från den tidigare ekvation (3). De andra 12 möjligheterna har en sannolikhet på

eller

. Då blir sannolikheten för den första nukleotiden

, då det existerar fyra möjligheter (A, T, C, G). Sanno- likheten att det sker en transition är p

, detta tillsammans blir då

. Är det istället en transversion blir det istället

och

. Det går enkelt att se hur ekvationerna (32) och (33) enbart skiljer sig på proportionerna av konstanter- na. Sannolikheten för maximum likelihood, L, är väldigt små och blir därför ganska besvärliga att arbeta med. Det är därför vanligt att man istället an- vänder logaritmen, `(d) = log {L(d)}. Logaritmfunktionen är monoton och den uppnår samma resultat som är L(d

) > L(d

) om och endast om `(d

) > `(d

) . Logaritmfunktionen blir följande

`(d; x) = log {L(d; x)} = x log( 1 16 − 1

16 e

)+(n−x)log( 1 16 + 1

16 e

) (34)

Genom att

= 0 kan man bestämma att logaritmen, `, är maximerad av maximum likelihood, L. Därifrån kan man få det skattade avståndet, b d ,

d = − ˆ 3

4 log(1 − 4 3

x

n ) (35)

Detta är precis lika som avståndet i ekvationen (7) för JC69-modellen som deriverades och hade i den ekvationen ˆp =

.

4.2 K80

K80 -modellen har era parametrar än JC69, vilket gör att maximum likelihood- metoden behöver justeras lite i jämfört med det tidigare avsnittet. Maximum likelihood-metoden tillämpas för att uppskatta sekvensens avstånd, d, samt transitions-/transversionsfrekvensens förhållande till, κ. Informationen som krävs för metoden är antalet nukleotider, n, samt antalet antalet övergångar, n

, och antalet transversionaler, n

. Sannolikheten beräknas för en konstant plats (ex- empelvis T T ) är

, och sannolikheten för en transitions skillnaden (exempel TC) är

. Till sist är sannolikheten för en transversell skillnad (till exempel TA)

. I ekvation (13) anges vad p

, p

och p

står för. Log-likelihood är

`(d, κ | n

, n

) = log{f (n

, n

| d, κ)}

=(n − n

− n

) log( p

4 ) + n

log( p

4 ) + n

log( p

4 )

Maximum likelihood-metoden av avståndet, d, och transitions- /transversions- frekvenskvoten, κ, härleds genom

= 0 och

= 0 . Detta kan lösas med ekvation (15), S =

och V =

.

5 Uppbyggnad av fylogenetiska träd

Det kommer att visas i en enklare grad hur man konstruerar och skapar fyloge- netiska träd. Ett fylogenetiska träd är en trädliknande graf där man ingående studerar relationer mellan arter, gener eller individer. Inom matematiken byggs grafer upp med hörn och kanter som bygger upp trädet bildligt. Här kommer det att skrivs om fylogeniträd uppbyggnad av arter. När ett träd konstrueras för arter kommer de externa noderna, eller löven som de också kallas, att re- presentera de arter som existerar idag och de interna noderna är arter som är utdöda. Slutligen vid roten av de fylogenetiska träd nns förfäderna. Träden ritas vanligtvis med roten högst upp där de interna noderna infaller nedanför och avslutas trädet med de externa noderna. Det är så generellt träden är upp- byggda, men det nns såklart era modeller utav uppbyggnad av träd. Ett träd kan till exempel vara orotad, då är man inte är säker på vilken förfadern är.

Det ger trädet en lite mer rundare form. Har man en evolutionsfrekvens som

är konstant över tiden så kallar man det för den molekylära klockan. Då av-

ståndsmatrisen och maximum likelihood-metoden identierar roten och sedan

konstruerar trädet. Utan den molekylära klockan är de svårt att identiera ro- ten, vilket gör det svårare att skapa ett träd. I ett fylogenetiskt träd brukar man kalla de närmsta besläktade arterna för ingrupper och de som är släkt på lite längre håll för utgrupper. Två vanliga trädtyper är kladogram och fylogram.

Ett kladogram är ett träd som inte visar någon information av grenarnas längd medan ett fylogram ger information genom grenarnas längder.

5.1 Avstånd mellan arter

Avståndsmetoder innebär två steg; beräkning av genetiska avstånd mellan två arter och rekonstruktion av ett fylogenetiskt träd från en avståndsmatris. En av de simplare avståndsmetoden är kanske UPGMA (Sokal och Sneath 1963).

Denna metod är baserad på den molekylära klockans antagande och genererar rotade träd. Det är tillämpligt på befolkningsuppgifter och används sällan för att analysera data arter, eftersom klockan ofta krävs när sekvenserna är avvi- kande. En metod som inte kräver den molekylära klockans antagande är den minstakvadratmetoden.

5.1.1 Minstakvadratmetoden

Minstakvadratmetoden tar den parvisa avståndsmatrisen, med given data, och uppskattar grenarnas längd på ett träd genom att para ihop deras avstånd så noga som möjligt. Det görs genom att minimera summan av kvadratskillnaden mellan det förutbestämda avståndet, d, och det skattade avståndet, ˆ d . De förut- bestämda avståndet beräknas genom att summera grenarnas längd mellan två gränsande arterna. Låt avståndet mellan art i och j vara d

. Låt summan av grenarnas längd från art i till j vara ˆ d

. Därefter används minstakvadratmeto- den för att minimera summan över samtliga par i och j med kvadratskillnaden (d

− ˆ d

)

, så att trädet passar avståndet så mycket som möjligt. Summan av kvadratskillnaden är som följer

S = X

i<j

(d

− ˆ d

)

(36)

Om man exempelvis har fyra arter a, b, c och d som bildar följande träd

Fig.4: Ett artträd som demonstrerar för minstakvadratmetoden där krite- rier för grenarna ((a, b), c, d)

Trädet består av fem grenar, t

, t

, t

, t

och t

. Det förutsagda avståndet mellan till exempel a och b är t

+ t

medan för a och c är avståndet t

+ t

+ t

. För a och d är avståndet t

+ t

+ t

, för b och c är avståndet t

+ t

+ t

och till sist är avståndet för c och d är lika med t

+ t

. Summan av kvadratskillnaden för detta fall är följande

S = X

i<j

(d

− ˆ d

)

= (d

− ˆ d

)

+(d

− ˆ d

)

+ (d

− ˆ d

)

+ (d

− ˆ d

)

+ (d

− ˆ d

)

+ (d

− ˆ d

)

5.2 Maximum likelihood-metoden - er generationer

Här diskuteras och beräknas sannolikheten för multipla sekvenser på ett fyloge- netiskt träd. Detta kommer att ske som en naturlig förlängning från de tidigare beräkningar av avståndet mellan två sekvenser. Boken går igenom två sätt att beräkna detta på, men i det här arbetet blir det bara fokus på en av möjliga metoder. Den metoden uppskattar parametrar i den evolutionära modellen och testar en hypotes om den evolutionära processen när en trädtopologi är känd samt xerad. Där tillämpas maximum likelihood-metoden som har många bra statistiska egenskaper och som ger en kraftfull och exibel för denna analys.

5.2.1 Likelihood beräkningar på träd

Som det tidigare förklarades denierar maximum likelihood-metoden sannolik-

heten för att observera data för en given parameter, även fast de anses vara en

funktion av parametrar. Här kommer man att utgå från K80 metoden. Man

utgår från att de olika platserna utvecklas oberoende av varandra och att en

gren är oberoende av en annan gren.

Fig.5: Ett träd med 5 arter som används för att demonstrera exemplet som används under maximum likelihood-funktionen. Grenarnas längd är mätta med förväntat antal nukleotid substitutioner per plats.

Förfäderna i trädet i Figur 5 är noderna 0, 6, 7 och 8, där 0 är roten. Grenarnas längder betecknas t

, där i står för noden grenen går till. Parametrarna i model- len inkluderar grenarnas längd och dess transition-/transversionsfrekvensen, κ , med gemensamma betäckningar θ = {t

, t

, t

, t

, t

, t

, t

, t

, κ} . Eftersom an- tagandet av en oberoende evolution mellan platserna där sannolikheten av hela datasekvensen är produkten av sannolikheten för enskild individs plats. Detta är ekvivalent med logaritmen av maximum likelihood-metoden är summan över platserna i sekvensen.

` = log(L) = X

log{f (x

| θ)}

6 Resultat

6.1 Människa (Homo sapiens) D38112 mot Schimpans (Pan troglodytes)

Här jämförs genfrekvensen D38112 hos en människa med en schimpans. Genfre-

kvensen nns under bilagor, bilaga 8.3.1. De först 960 nukleotiderna har blivit

jämförda med varandra, dock har de nukleotiderna som inte kunnat jämföras

tagits bort så det totala antalet jämförda nukleotider är 956 stycken. Tabell 1

är lätt avrundad med små modieringar.

Människa

Schimpans T C A G P π

T 195(0.2039) 15(0.0157) 0(0) 4(0.0042) 0.2238 C 9(0.0094) 247(0.2584) 2(0.0021) 0(0) 0.2699

A 0(0) 0(0) 309(0.3232) 0(0) 0.3232

G 0(0) 0(0) 6(0.0063) 169(0.1768) 0.1831

P π 0.2133 0.2741 0.3316 0.1810 P 1

Tabell 1.

Tabell 1 visar alla nukleotider i sekvensen och vilka nukleotider som skiljer sig åt mellan arterna. I parenteserna visas dierensen för alla möjliga utfall. Det genomsnittliga värdet på frekvensen för de olika nukleotiderna är T = 0.21855, C = 0.24685 , A = 0.30075 och G = 0.18205.

6.1.1 JC69 modellen

De totala antalet nukleotider från genfrekvensen som är jämförda i det här fallet är n = 956. De nukleotider som skiljer sig från varandra summeras ihop, x = 2+

6+9+15+4 = 36 . Tillsammans med n och x räknas andelen olika platser ut, ˆp =

x

n

=

= 0.03765690377 . De skattade avståndet beräknas från ekvation (7), d = − ˆ

log(1−

p) = 0.03863515 ˆ . Därefter räknas variansen ut från ekvation (8), var( ˆ d) =

3p)ˆ²

= 0.00004202 . Variansen behövs för att kunna räkna ut standardfelet, som är roten ur på variansen, så standardfelet blir ε = q

var( ˆ d) = 0.0064823 . Till sist tillämpas ett approximerat 95%-kondensintervall

d ± λ ˆ

· ε =

( 0.0513405 0.0259298

där λ

= 1.96 , värdet nns i tabellen under bilagor, bilaga 8.1. Om variansen istället räknas ut med sannolikheten, p, från ekvation (8), var(ˆp) =

= 0.0000379 . Det nya standardfelet räknas ut till ε = pvar(ˆp) = 0.006157 och de approximerade 95%-kondensintervallet blir

d ± λ ˆ

· ε =

( 0.0507026 0.0265677

Från ekvation (33) kan man räkna ut maximum liklihood där p : ˆp =

`( ˆ d) =`(ˆ p) = x log( x

12n ) + (n − x) log( n − x 4n )

=36 · log( 36

12 · 956 ) + (956 − 36) log( 956 − 36

4 · 956 ) = −1518.213558

Genom att sänka värdet av log-likelihood till χ

från toppen skapas ett approximerat 95%-kondensintervall för avståndet, d, och sannolikheten, p. Den generella formen är χ

där k står för grad av frihet och det står för antalet parametrar. Då får vi χ

/2 = 3.841/2 = 1.921 , vilket ger `(p) = `(d) = `( ˆ d)−

χ

/2 = −1520.134559 . Därefter görs en fplot av maximum likelihood- funktionen av avståndet, d. Därefter studeras intervallet för avståndet, d, under χ

som blir `(d) = −1520.134559.

Fig 6: Graf av loglikelihood-funktionen över avståndet, d, för JC69. Grafen visar vilka värden kondensintervallet under χ

antar.

Grafen ger oss ett intervall för avståndet, d, på (0.02659, 0.05104). Gör man en fplot på maximum likelihood av sannolikheten, p, istället och studerar inter- vallet för det under χ

, `(p) = −1520.134599.

Fig 7: Graf av loglikelihood-funktionen över sannolikheten, p, för JC69.

Grafen visar vilka värden kondensintervallet under χ

antar.

Intervallet över sannolikheten, p, blir som följer (0.02714, 0.0506).

P (t) = e

=









p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t)







 där övergångssannolikheterna blir

( p

(t) =

+

e

=

+

e

= 0.962343 p

(t) =

−

e

−

e

= 0.0125523 En kontroll görs, p

+ 3 · p

= 1 , och det visar att det stämmer.

6.1.2 K80-modellen

Här beräknas K80-modellen för samma data som för JC69. Beräkningarna för K80 blir bara lite mer precisare än för JC69. Istället för att endast räkna alla nukleotider som skiljer sig från varandra delar man upp de i två grup- per, transitions och transversions skillnad. Andelen av alla transitions skill- nader blir S = (9 + 15 + 6 + 0)/956 = 30/956 = 15/478 och transversions skillnaden blir V = (0 + 0 + 0 + 0 + 2 + 4 + 0 + 0 + 0)/956 = 6/956 = 3/478 . Efter detta beräknas de skattade avståndet ut från ekvation (14), ˆ d =

−

log(1 − 2S − V ) −

log(1 − 2V ) = 0.03892616 . Variansen av det skatta- de avståndet är var( ˆ d) =

n

, där a = (1 − 2S − V )

och b =

(1 − 2S − V )

+ (1 − 2V )



. Genfrekvensens data räknas ut och blir a = 1.074157 , b = 1.043435 och var( ˆ d) = 0.0000434067 . Med all insamlad data kan det approximerade 95%-kondensintervallet beräknas till

d ± λ ˆ

· ε =

( 0.05183938 0.02601294 där standardelet är ε = q

var( ˆ d) = 0.00658838 . Från ekvation (15) kan den skattade transitions- och transversionsfrekvenskvoten beräknas, som blir ˆκ =

2log(1−2S−V )

log(1−2V )

− 1 = 10.326454 . Detta visar oss att transitionsfrekvenskvoten är ungefär 10 gånger högre än transversionfrekvenen.

c αt = − 1

2 log(1 − 2S − V ) + 1

4 log(1 − 2V )

= − 1

2 log(0.9309623) + 1

4 log(0.9874477) = 0.0326103

d 2βt = − 1

2 log(1 − 2V ) = − 1

2 log(0.9874476) = 0.00631587

Maximum likelihood-metoden beräknas enligt ekvation (36) och ser ut som följer

`(d, κ | n

, n

) = (n − n

− n

) log( p

4 ) + n

log( p

4 ) + n

log( p

4 ) där övergångssannolikheten blir



 

 

p

(t) =

+

e

+

e

=

+

e

+

e

= 0.962343 p

(t) =

+

e

−

e

=

+

e

−

e

= 0.031381

p

(t) =

−

e

=

−

e

= 0.003138 Detta ger oss följande

`(d, κ | n

, n

) = (956−30−6) log( p

4 )+30 log( p

4 )+6 log( p

4 ) = − 1468.654387 En kontroll kan görs för att kotrollera att värdena för övergångssannolikheter- na stämmer, p

+ p

+ 2 · p

= 1 . I vårt fall stämmer kontrollera. Även här studeras χ

med en frihetsgrad 1 och från toppen skapas ett approxime- rat 95%-kondensintervall för avståndet, d. Där χ

/2 = 1.921 vilket ger oss

`( ˆ d) − χ

/2 = −1500.963602 .

Fig 8: Graf av likelihood över sannolikheten, p, för K80. Grafen visar vilka värden kondensintervallet under χ

antar.

Från grafen får vi intervall (0.02903, 0.05089). För att kunna räkna ut variansen av ˆ d över ˆκ behövs bland annat variansen av S över V samt jacobianen, J.

Ekvation (19) visar hur formeln ser ut. Variansen av S över V beräknas på följande sett

var( S V ) =

S(1−S)

n

−

−

!

=

 0.00003179498 −0.00000020601

−0.00000020601 0.00000652381



Där jacobianen räknas ut såhär

J =

1 1−2S−V

1

2(1−2V )

+

−4

(1−2S−V ) log(1−2V )

−2

(1−2S−V )log(1−2V )

+

!

=

 1.0741573 1.0434346 340.14539 −1646.05531



Med all denna information kan nu variansen av det skattade avståndet, ˆ d , över det skattade transitions-/transversionsfrekvenskvoten,ˆκ, beräknas

var( d ˆ ˆ

κ ) = J · var( S

V ) · J

=

 0.0000433265 0.0007030852 0.0007030852 21.585582163



6.2 Människa (Homo sapiens) D38112 mot Gorilla (Goril- la gorilla)

Här jämförs genfrekvensen D38112 hos en människa med en gorilla. De först 960 nukleotiderna har blivit jämförda med varandra, dock har de nukleotiderna som inte kunnat jämföras tagits bort så det totala antalet jämförda nukleotider är 956 nukleotider. Genfrekvenserna nns under bilagor, bilaga 8.3.2.

Människa

Gorilla T C A G P π

T 194(0.2029) 12(0.0126) 0(0) 6(0.0063) 0.2218 C 12(0.0126) 244(0.2552) 1(0.0010) 0(0) 0.2688

A 0(0) 3(0.0031) 312(0.3264) 0(0) 0.3295

G 0(0) 0(0) 5(0.0052) 167(0.1747) 0.1799

P π 0.2155 0.2709 0.3326 0.181 P 1

Tabell 2.

Tabell 2, visar alla nukleotider i sekvensen och vilka nukleotider som skiljer sig åt mellan arterna. I parenteserna visas dierensen för alla möjliga utfall.

Medelvärdet av frekvensen för de olika nukleotiderna är T = 0.21865, C = 0.26985 , A = 0.33105 och G = 0.18045.

6.2.1 JC69 modellen

De totala antalet nukleotider från genfrekvensen som är jämförda är n = 956.

De nukleotider som skiljer sig från varandra summeras ihop x = 1 + 5 + 12 + 3 + 12 + 6 = 39 . Tillsammans med n och x räknas andelen olika platser ut, ˆp =

=

39

956

= 0.04079498 . Det uppskattade avståndet beräknas genom ekvation (7),

d = − ˆ

log(1 −

p) = 0.041946 ˆ . Därefter räknas variansen ut för det skattade

avståndet, ˆ d , från ekvation (8) vilket ger oss

3p)ˆ²

= 0.00004578 . Vari- ansen behövs för att kunna beräkna standardfelet, detta görs genom att roten ur på variansen, ε = q

var( ˆ d) = 0.0067658 . Till sist tillämpas ett approximerat 95%-kondensintervall

d ± λ ˆ

· ε =

( 0.0552074 0.0286854

där λ

= 1.96 och tabellen nns under bilagor, bilaga 8.1. Om variansen kalkyleras ut med sannolikheten, p, istället i ekvation (8) får vi var(ˆp) =

ˆ p(1− ˆp)

n

0.00004093 . Det nya standardfelet räknas ut och blir ε = pvar(ˆp) = 0.0063978 och det approximerade 95%-kondensintervallet justeras till följande

d ± λ ˆ

· ε =

( 0.0544860 0.0294067 Maximum likelihood ger oss

`( ˆ d) =`(ˆ p) = x log( x

12n ) + (n − x) log( n − x 4n )

=39 · log( 39

12 · 956 ) + (956 − 39) log( 956 − 39

4 · 956 ) = −1531.105397

Genom att sänka log-likelihood till χ

från maximum likelihood-metoden skapas ett 95% intervall för avståndet, d, och sannolikheten, p. Där χ

/2 = 3.841/2 = 1.921 vilket ger `(p) = `(d) = `( ˆ d) − χ

/2 = −1533.026397 . Där- efter görs en fplot av maximum likelihood funktionen av avståndet,d. Därefter studeras intervallet av avståndet, d, under χ

alltså för `(d) = −1533.026397.

Fig 9: Graf av likelihood över avståndet, d, för JC69. Grafen visar vilka

värden kondensintervallet under χ

antar.

Detta ger oss ett intervall för avståndet, d, på (0.02956, 0.05715). Om en fplot istället görs på maximum likelihood av sannolikheten, p, och studerar intervallet av p under χ

, `(p) = −1533.026397.

Fig 10: Graf av likelihood över sannolikheten, p, för JC69. Grafen visar vilka värden kondensintervallet under χ

antar.

Intervallet över sannolikheten, p, blir som följer (0.03424, 0.05337).

P (t) = e

=









p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t) p

(t)







 där övergångssannolikheten blir

( p

(t) =

+

e

=

+

e

= 0.9592050 p

(t) =

−

e

=

−

e

= 0.01359832 En kontroll görs, p

+ 3 · p

= 1 , och det stämmer.

6.2.2 K80-modellen

Här beräknas K80-modellen istället med samma data som innan. Istället för att beräkna alla nukleotider som skiljer sig från varandra delar man upp dem i två grupper, transitions och transversions skillnader. Andelen av transitions skillnaden blir S = (12 + 12 + 5 + 0)/956 = 29/956

och transversions skillnaden blir V = (0 + 0 + 3 + 0 + 1 + 6 + 0 + 0)/956 = 10/956 = 5/478 . Efter detta kalkyleras det skattade avståndet ut från ek- vation (14) ˆ d = −

log(1 − 2S − V ) −

log(1 − 2V ) = −

log(0.928870) −

1

4

log(0.979079) = 0.0421787 . Variansen av det skattade avståndet är var( ˆ d) =

a²S+b²V −(aS+bV )²

n

, där a = (1−2S−V )

och b =

(1 − 2S − V )

+ (1 − 2V )

 .

Med genfrekvensens data blir a = 1.076577, b = 1.048972 och

var( ˆ d) = 0.00004690 . Med all insamlad data kan ett approximerat 95%- kondensintervall beräknas till

d ± λ ˆ

· ε =

( 0.05560158 0.02875581 där standardfelet är ε = q

var( ˆ d) = 0.00684841 . Från ekvation (15) kan den skattade transitions- och transversionsfrekvenskvoten beräknas till ˆκ =

log(1−2V )

− 1 = 5.9799117 . Detta visar oss att transitionsfrekvensen är nästan 6 gånger hög- re än transversionfrekvensen.

c αt = − 1

2 log(1−2S−V )+ 1

4 log(1−2V ) = − 1

2 log(0.928870)+ 1

4 log(0.979079) = 0.0316075

d 2βt = − 1

2 log(1 − 2V ) = − 1

2 log(0.979079) = 0.0105712 Maximum likelihood-metoden beräknas enligt ekvation (35)

`(d, κ | n

, n

) = (n − n

− n

) log( p

4 ) + n

log( p

4 ) + n

log( p

4 ) där övergångssannolikheten blir



 

 

p

(t) =

+

e

+

e

=

+

e

+

e

= 0.95920502 p

(t) =

+

e

−

e

=

+

e

−

e

= 0.03033473

p

(t) =

−

e

=

−

e

= 0.0052301 Vilket ger oss

`(d, κ | n

, n

) = (956−29−10) log( p

4 )+29 log( p

4 )+10 log( p

4 ) = −1471.895889

Kontrollerar övergångssannolikheternas värden som innan, p

+ p

+ 2 · p

= 1 ,

och det stämmer. Studera χ

med frihetsgraden 1 och från maximum skapas

ett approximerat 95%-kondensintervall för avståndet, d. Där χ

/2 = 1.921

ger `( ˆ d) − χ

/2 = −1519.313464 .

Fig 11: Graf av likelihood över avståndet, d, för K80. Grafen visar vilka värden kondensintervallet under χ

antar.

Intervallet för avståndet, d, är (0.03086, 0.0545). För att kunna räkna ut varian-

sen av det skattade avståndet, ˆ d, över det skattade transitions-/transversionsfrekvenskvoten, ˆ

κ, behövs bland annat variansen av S över V samt jacobianen, J. Ekvation (19) visar hur formeln ser ut. Variansen av S över V beräknas på här

var( S V ) =

S(1−S)

n

−

−

!

=

 0.000030768 −0.000000332

−0.000000332 0.0000108272



Där jacobianen beräknas såhär

J =

1 1−2S−V

1

2(1−2V )

+

−4 (1−2S−V )log(1−2V )

−2

(1−2S−V )log(1−2V )

+

!

=

 1.076577 1.04897 203.6807 −572.5431



Med all denna information kan nu variansen av det skattade avståndet, ˆ d , över det skattade transitions-/transversionsfrekvenskvoten, ˆκ, beräknas till

var( d ˆ ˆ

κ ) = J · var( S

V ) · J

=

 0.00004683 0.00037784 0.00037785 4.90309071



6.3 Människa (Homo sapiens) D38112 mot Bonobo (Pan paniscus)

Här jämförs genfrekvensen D38112 hos en människa med en schimpans. De först

960 nukleotiderna har blivit jämförda med varandra, dock har de nukleotiderna