BAB III
3.1.
KOMBINASI MMT DENGAN FILOGENETIK
Pendahuluan
Aplikasi MMT dalam bioinformatik pada dasarnya sangat banyak. Salah satunya adalah dalam filogenetik. Filogenetik adalah cabang ilmu yang mempelajari evolusi biologi. Pada umumnya, model-model dalam filogenetik memiliki rate evolusi yang tetap untuk setiap site, padahal dalam kenyataannya tidak seperti itu. Oleh karena itu dalam tesis ini dibangun model filogenetik dengan rate evolusi yang berbeda untuk setiap site. Karena rate evolusi di setiap site bersifat tersembunyi, maka dilibatkan model MMT. Dalam tesis ini diasumsikan rate evolusi mempunyai distribusi Gamma dengan parameter α dan β merupakan fungsi dari intensitas perubahan basanya. Kategori rate ditentukan dari sampel acak yang dibangkitkan. Dalam tesis ini dianggap ada 3 kategori rate evolusi, sehingga terdapat 3 keadaan yang tersembunyi. Pada MMT, peluang emisi bisa ditentukan dari distribusi meltinomial sebagai proporsi masing-masing elemen observasi untuk setiap keadaan. Dalam banyak referensi biasanya proporsinya sudah diketahui. Yang menjadi pertanyaan adalah: ”informasi proporsi tersebut didapat dari mana?”padahal keadaannya tersembunyi. Oleh
karena
itu
banyak
peneliti
yang
menentukan
peluang
emisi
dari
informasi/penelitian sebelumnya. Ada juga yang menentukan peluang emisinya secara acak. Alasannya, meskipun dipilih secara acak, arah kekonvergenannya akan tetap sama, yang beda adalah laju kekonvergenannya. Tetapi alasan ini belum bisa dibuktikan secara pasti, jadi kita tidak bisa percaya sepenuhnya. Peluang emisi dalam tesis ini ditentukan dari kategori rate evolusi yang dihasilkan dari rate evolusi yang berdistribusi Gamma. Sedangkan matriks transisi di sini menyatakan peluang dari kombinasi rate evolusi yang berbeda untuk setiap site dan melibatkan fungsi otokorelasi dari rate evolusi untuk site yang berdekatan. Oleh karena itu, sebelum kita mempelajari lebih jauh tentang kombinasi model filogenetik dengan MMT, harus dipelajari terlebih dahulu bagaimana membangun model filogenetik.
35
3.2.
Model Filogenetik
Secara umum, model filogenetik didefinisikan sebagaiψ (τ , π , Q, b) , sebagai 4tuple yang terdiri atas topologi pohon τ , frekuensi basa π , matriks intensitas substitusi Q, dan panjang cabang b. Agar lebih jelas lagi, ditampilkan flow-chart dari kerangka filogenetik seperti yang tampak di bawah ini:
Memilih Barisan
Butuh Perbaikan
Mengatur Barisan
Model Evolusi
Memilih Model dan Metode; Membangun Pohon
Metode
Membangun Pohon
Evaluasi Pohon Good
Interpretasi Filogenetik
Berdasarkan flow-chart di atas, setelah memilih dan mengatur barisan DNA, tahapan selanjutnya adalah pemilihan model. Dalam pemilihan model ψ (τ , π , Q, b) , berarti kita harus menjelaskan satu-satu parameter yang mempengaruhi model tersebut. Ada 4 metode yang sering digunakan, diantaranya adalah metode parsimoni, distance, Bayes, dan maksimum likelihood. Dalam tesis ini yang akan digunakan adalah metode maksimum likelihood. Alasan pemilihan metode ini adalah karena estimator yang dihasilkan sudah dijamin bersifat konsisten. Beda misalkan dengan metode parsimoni yang memiliki estimator yang tidak konsisten untuk pohon evolusinya (Felsenstein, 1978). Dalam prakteknya, metode maksimum likelihood dan parsimoni memiliki hasil yang sama juga, dan jika dibuat peluang awalnya sama, maka metode maksimum laikelihood tidak lain adalah parsimoni juga. Yang disebut konsisten dalam hal ini adalah meminimumkan variansi. Hubungannya dengan topologi pohon, disebut konsisten jika ada beberapa orang yang mengambil sampel DNA yang sama, hasil pohon dan panjangnya tidak akan jauh berbeda. Untuk mengetahui topologi pohon, banyak peneliti yang menggunakan
36
metode resampling seperti bootstrap dan jacknife, dan hasilnya akan konsisten jika setelah dilakukan resampling, topologi pohonnya tetap saja seperti itu. Setelah model filogenetik terbangun, model tersebut kemudian dievaluasi. Proses evaluasi dilakukan dengan pendekatan goodness of fit. Namun sering kali para peneliti mengevaluasi model ini dengan menggunakan Akaike Information Criterion (AIC) atau Bayes Information Criterion (BIC). Prinsip AIC dan BIC hampir sama dengan meminimumkan variansi, artinya parameter-parameter yang dihasilkan dari proses estimasi akan dipilih untuk yang memiliki variansi yang minimum. Jika model sudah baik maka model tersebut bisa langsung digunakan, jika masih kurang baik, maka perlu dilakukan identifikasi model lagiseperti yang sebelumnya Namun dalam tesis ini tidak akan dibahas tentang evaluasi pohon, karena fokus kita ke MMT-nya. 3.3.
Topologi Pohon ( τ )
Topologi pohon berkaitan erat dengan bentuk pohon filogenetik. Pohon filogenetik terdiri atas titik cabang yang masing-masing titiknya menyatakan spesies dan cabangnya menyatakan hubungan antar spesies-spesies tersebut. Setiap cabang hanya menghubungkan dua buah titik yang berdekatan. Di bawah ini digambarkan pohon filogenetik beserta dengan komponen yang mempengaruhinya, Spesies puncak (ancestral)
Cabang Titik cabang
Panjang cabang
Spesies 1
Spesies 2
Spesies 3
Spesies 4
Spesies 5
Gambar 3 Komponen Pohon Filogenetik
Pohon filogenetik terdiri atas: titik cabang, cabang, puncak dan daun. Pohon filogenetik ini dibangun berdasarkan spesies yang ada pada daun. Hanya barisan DNA dari spesies-spesies ini yang diketahui sedangkan informasi DNA pada spesies di titik cabang dan puncak pohon tidak diketahui. Panjang cabang merepresentasikan durasi atau lamanya evolusi dari spesies tersebut.
37
Terdapat dua jenis pohon filogenetik yang sering digunakan, yaitu pohon berpuncak dan yang tidak berpuncak. Perbedaan dari kedua topologi tersebut adalah ada tidaknya spesies puncak. Dari kedua topologi pohon filogenetik tersebut, masingmasing jenis pohon dikelompokkan lagi ke dalam dua kelompok berdasarkan panjang cabang yaitu pohon berskala dan tidak berskala. Arti dari topologi pohon yang tidak berpuncak dengan cabang berskala adalah bahwa cabang-cabang pada pohon filogenetik tersebut dibuat berdasarkan skala sehingga panjang masing-masing cabang berbeda-beda. Satu satuan panjang cabang bisa berupa waktu (tahun, bulan dan sebagainya) atau berupa banyaknya perubahan basa yang terjadi. Sedangkan pada topologi pohon dengan panjang cabang tidak berskala, panjang cabangnya dibuat sama. Meski panjang cabang di gambar dengan ukuran yang sama, tidak berarti durasi atau waktu evolusi antar spesies adalah sama. Biasanya untuk pohon dengan panjang cabang tidak berskala, panjang cabang ditulis di atas cabangnya sebagai informasi. Berikutnya akan dibahas kombinasi topologi pohon filogenetik. Ada berbagai kemungkinan pohon filogenetik yang dapat dibangun. Banyaknya pohon filogenetik ini berdasarkan banyaknya spesies pada ujungnya. Sebagai contoh untuk banyak adalah spesies 4, maka terdapat 5 kombinasi pohon filogenetik berpuncak yang mungkin. Ada juga topologi pohon binari, yaitu topologi pohon dimana setiap titik cabangnya hanya membawahi 2 titik cabang lainnya dan jika diberikan N spesies maka banyaknya titik dalam pohon adalah 2N-1 dan banyaknya cabang adalah 2N-2. 3.4.
Frekuensi Basa ( π ) dan Matriks Intensitas Substitusi (Q)
Proses stokastik merupakan suatu proses dimana proses yang akan datang atau selanjutnya hanya tergantung pada proses yang sekarang dan tidak dipengaruhi oleh proses yang lalu. Dalam proses stokastik ada kaitan antara rantai Markov dan proses Poisson khususnya dalam kelahiran dan kematian murni. Proses perubahan basa nukleotid juga merupakan proses stokastik karena perubahan basa nukleotid pada spesies saat ini tergantung pada keadaan basa spesies anchestralnya yang hidup s sebelumnya. Semakin besar perubahan basa, semakin cepat perubahan spesiesnya. Dan semakin panjang, mutasi akan semakin cepat. Masa evolusi suatu spesies dari spesies anchestralnya merupakan masa yang sangat lama dan dalam waktu yang lama tersebut kita tidak dapat mengetahui berapa 38
banyak perubahan basa nukleotid yang terjadi. Yang ada hanyalah informasi mengenai basa nukleotid spesies tersebut tetapi tidak untuk basa nukleotid spesies anchestralnya. Peluang dalam selang waktu yang pendek terjadi lebih dari 1 perubahan sangat kecil dan dapat diabaikan. Oleh karena itu wajar jika banyaknya perubahan di sini memiliki distribusi Poisson dan proses perubahan tersebut merupakan proses Poisson. Misalkan M adalah banyaknya perubahan basa yang terjadi selama waktu s, s ≥ 0 dan mengikuti proses Poisson, dengan peluang didefinisikan sebagai berikut P(M = m ) = exp(− θs )
(θs )m , m = 0,1,2,K m!
(3.1)
dengan θ > 0 adalah banyaknya kejadian dalam satuan waktu/lokasi. Dalam kasus ini, θ adalah banyaknya perubahan basa per satuan waktu (disebut juga intensitas substitusi). Biasanya satuan dari θ adalah mya (million years ago). Jika S menyatakan waktu antar perubahan basa pada spesies u dari spesies anchestralnya, maka S berdistribusi eksponensial dengan parameter θ , P(S = s ) = exp(− θs ) , s ≥ 0
(3.2)
Bukti Misal: {M(s), s ≥ 0} suatu proses Poisson yang menyatakan banyaknya perubahan basa selama waktu s dan Si menyatakan waktu antara kejadian ke-(i – 1) dan ke–i, maka P(S1 > s ) = P(M (s ) = 0) = exp(− θs ) dan P(S 2 | S1 = t ) = P(tidak ada kejadian dalam selang (t , t + s ] | S1 = s ) = P(tidak ada kejadian dalam selang (t , t + s ]) = exp(− θs ) dan demikian seterusnya, dapat disimpulkan bahwa Si, i = 1, 2, … identik dan saling bebas memiliki distribusi eksponensial dengan parameter θ > 0 (Ross, 1996). Misal terdapat m perubahan dari spesies anchestral ke spesies u selama s waktu, maka peluang transisi m langkah dari basa i menjadi basa j sama dengan peluang transisi dari basa i menjadi basa j selama s waktu. Di sini ada dua hal
39
berbeda. Di satu sisi, kita berbicara dengan indeks parameter diskrit, namun di sisi lain kita berbicara dengan indeks parameter waktu (kontinu). Misal peluang pada saat s berada di basa j jika pada saat 0 berada di basa i atau peluang perubahan dari basa i ke basa j selama s adalah bij (s ) dan M adalah peubah acak yang menyatakan banyaknya perubahan basa dari basa i ke basa j berdistribusi Poisson dengan parameter θs , maka
[ ]
bij (s ) = E rij(M ) ∞
= ∑ rij(m ) P(M = m ) m =0 ∞
(θs )m
m =0
m!
= ∑ rij(m ) exp(− θs ) ∞
= exp(− θs )∑ rij
(m )
(θs )m m!
m=0
untuk i, j ∈ {A, C , G, T }
Atau jika ditulis dalam bentuk matriks,
⎛ bAA (s ) ⎜ ⎜ bCA (s ) ⎜ b (s ) ⎜ GA ⎜ b (s ) ⎝ TA
bAC (s ) bAG (s ) bCC (s ) bCG (s )
bGC (s ) bGG (s ) bTC (s ) bTG (s )
m ∞ ⎛ ( m ) (θs ) ⎜ exp(− θs )∑ rAA m! ⎜ m =0 m ∞ bAT (s )⎞ ⎜ ( θ m ( ) ⎟ ⎜ exp(− θs )∑ rCA s ) bCT (s ) ⎟ ⎜ m! m =0 = m ∞ bGT (s )⎟ ⎜ ( m ) (θs ) ⎟ ⎜ exp(− θs )∑ rGA ⎟ m! bTT (s ) ⎠ ⎜ m =0 m ∞ ⎜ ( m ) (θs ) ⎜ exp(− θs )∑ rTA m! m =0 ⎝
⎛ ∞ (m ) (θs )m ⎜ ∑ rAA m! ⎜ m =0 ⎜ ∞ (m ) (θs )m ⎜ ∑ rCA m! = exp(− θs )⎜ m∞=0 m ⎜ ( m ) (θs ) ⎜ ∑ rGA m! ⎜ m =0 m ∞ ⎜ ( m ) (θs ) ⎜ ∑ rTA m! ⎝ m =0
∞
(m ) exp(− θs )∑ rAC m =0 ∞
exp(− θs )∑ rCC
(m )
(θs )m m!
(θs )m
m! m (m ) (θs ) exp(− θs )∑ rGC m! m =0 m ∞ (m ) (θs ) exp(− θs )∑ rTC ! m m =0 m =0 ∞
∞
∑ rAC(m )
(θs )m
m! m ( m ) (θs ) r ∑ CC m! m =0 m ∞ ( m ) (θs ) ∑ rGC m! m =0 m ∞ ( m ) (θs ) ∑ rTC m! m =0
m =0 ∞
∞
(m ) exp(− θs )∑ rAG m=0 ∞
exp(− θs )∑ rCG
(m )
(θs )m
∞
(m ) exp(− θs )∑ rAT
m!
m=0 ∞
(θs )m
exp(− θs )∑ rCT
m! m (m ) (θs ) exp(− θs )∑ rGG m! m=0 m ∞ (m ) (θs ) exp(− θs )∑ rTG ! m m=0 m=0 ∞
∞
∑ rAG(m )
(θs )m
m! m ( m ) (θs ) r ∑ CG m! m =0 m ∞ ( m ) (θs ) ∑ rGG m! m =0 m ∞ ( m ) (θs ) ∑ rTG m! m =0
m =0 ∞
m=0 ∞
(m )
exp(− θs )∑ rGT
(θs )m ⎞⎟ m! ⎟
(θs )m ⎟⎟ m! ⎟
(m )
(θs )m ⎟
⎟ m! ⎟ m (θs ) ⎟ exp(− θs )∑ rTT(m ) m! ⎟⎠ m=0 m=0 ∞
∞
∑ rAT(m )
(θs )m ⎞⎟
m! m ( m ) (θs ) r ∑ CT m! m= 0 m ∞ ( m ) (θs ) rGT ∑ m! m= 0 m ∞ ( m ) (θs ) rTT ∑ m! m= 0 m=0 ∞
⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
40
⎛ (m ) (θs )m ⎜ rAA m! ⎜ ⎜ (m ) (θs )m ∞ ⎜ rCA m! = exp(− θs )∑ ⎜ m m = 0 ⎜ ( m ) (θs ) r ⎜ GA m! ⎜ (θs )m ⎜ rTA(m ) m! ⎝
m! m ( m ) (θs ) rGC m! m ( θ (m ) s ) rTC m!
m! m ( m ) (θs ) rGG m! m ( θ (m ) s ) rTG m!
(m ) ⎛ rAA ⎜ (m ) ∞ θs )m ⎜ rCA ( = exp(− θs )∑ ⎜ r (m ) m = 0 m! ⎜ GA ⎜ r (m ) ⎝ TA
(m ) rAC (m ) rCC (m ) rGC (m ) rTC
(m ) ⎞ rAT ⎟ (m ) rCT ⎟ (m ) ⎟ rGT ⎟ rTT(m ) ⎟⎠
(m ) rAC
(θs )m m!
m ( m ) (θs ) r CC
(m ) rAG (m ) rCG (m ) rGG (m ) rTG
(m ) rAG
(θs )m m!
m ( m ) (θs ) r CG
(m ) rAT
(θs )m ⎞⎟ m! ⎟
m ( m ) (θs ) ⎟ r CT
m! m ( m ) (θs ) rGT m! m ( s) θ rTT(m ) m!
⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
θs )m m ( = exp(− θs )∑ R ∞
m!
m =0 ∞
(Rθs )m
m=0
m!
= exp(− θs )∑
= exp(− θs ) exp(Rθs ) = exp(Rθs − θs ) = exp((R − I )θs )
Jadi
B(s ) = exp((R − I )θs ) , θ > 0
(3.3)
Perhatikan bahwa exp((R − I )θs ) bukan berarti setiap elemen dari matriks (R-
I)θs ditransformasikan dengan fungsi eksponensial. Eksponensial dari matriks D didefinisikan sebagai
Dk k =0 k! ∞
exp(D ) = ∑
Sekarang kita telah memiliki hubungan antara Rantai Markov diskrit dengan matriks peluang transisi R dan Rantai Markov kontinu dengan matriks peluang transisi B(s) sebagai berikut
41
B(s ) = exp((R − I )θs ) = exp(Qs ) Matriks B(s) di atas merupakan matriks peluang transisi dari Rantai Markov kontinu. Jadi sekarang kita dapat menggunakan rantai Markov kontinu dalam melihat perubahan basa sepanjang s satuan waktu. Jadi sekarang dalam kasus filogenetik, kita melihat perubahan basa dari spesies anchestral ke spesies u sepanjang s satuan waktu sebagai Rantai Markov kontinu dengan peluang transisi P(s) tidak lagi perubahan basa hanya dari spesies ancestral ke spesies saat ini saja. Perhatikan kembali matriks (R − I )θ pada persamaan (3.3) di atas, ⎧⎛ rAA ⎪⎜ r (R − I )θ = ⎪⎨⎜⎜ CA ⎪⎜ rGA ⎪⎜⎝ rTA ⎩
rAC
rAG
rCC
rCG
rGC
rGG
rTC
rTG
rAT ⎞ ⎛ 1 ⎟ ⎜ rCT ⎟ ⎜ 0 − rGT ⎟ ⎜ 0 ⎟ ⎜ rTT ⎟⎠ ⎜⎝ 0
0 0 0 ⎞⎫ ⎟⎪ 1 0 0 ⎟⎪ ⎬θ 0 1 0 ⎟⎪ ⎟ 0 0 1 ⎟⎠⎪⎭
rAG rAT ⎞ ⎛ rAA − 1 rAC ⎜ ⎟ rCC − 1 rCG rCT ⎟ ⎜ rCA θ =⎜ rGA rGC rGG − 1 rGT ⎟ ⎜ ⎟ ⎜ r rTC rTG rTT − 1⎟⎠ ⎝ TA ⎛ − (rAC ⎜ ⎜ =⎜ ⎜ ⎜ ⎝
+ rAG + rAT ) − (rCA rCA rGA rTA
rAC + rCG + rCT ) − (rGA rGC rTC
rAG rCG + rGC + rGT ) − (rTA rTG
rAT ⎞ ⎟ rCT ⎟ ⎟θ rGT ⎟ + rTC + rTG )⎟⎠
Misalkan ⎛ − (rAC ⎜ ⎜ Q=⎜ ⎜ ⎜ ⎝
+ rAG + rAT ) − (rCA rCA rGA rTA
rAC + rCG + rCT ) − (rGA rGC rTC
rAG rCG + rGC + rGT ) − (rTA rTG
rAT ⎞ ⎟ rCT ⎟ ⎟θ rGT ⎟ + rTC + rTG )⎟⎠
maka qii = −θ ∑ rij adalah banyaknya perubahan basa dari basa i per satuan waktu j j ≠i
atau intensitas perubahan basa (intensitas substitusi). Matriks Q ini yang kemudian disebut sebagai matriks intensitas substitusi.
42
Terdapat beberapa macam model matriks intensitas substitusi yang sering digunakan dalam kasus filogenetik, antara lain: a. Model Jukes dan Cantor Model ini diperkenalkan pertama kali oleh Jukes dan Cantor pada tahun 1969. Model ini merupakan model yang paling sederhana yaitu hanya menggunakan satu parameter. Asumsi yang digunakan di dalam model ini adalah substitusi antara keempat basa nukleotid terjadi secara acak. Dengan kata lain peluang terpilihnya masing-masing basa nukleotid seragam. Dalam model ini, intensitas substitusi untuk masing-masing transisi yang mungkin adalah α , konstan. Perhatikan diagram transisi berikut ini.
Frekuensi basa pada saat setimbang untuk model ini adalah sama
(π A = π C = π G = π T = 14 ) dan matriks intensitas substitusi dari model ini
didefinisikan sebagai berikut, ⎛ − 34α ⎜ α ⎜ Q = ⎜ α4 ⎜ 4 ⎜ α ⎝ 4
α
−
4 3α 4
α
4
α
4
α
4
− 34α
α
α
4
4
α
⎞ ⎟ α 4 ⎟ α ⎟ 4 ⎟ − 34α ⎟⎠ 4
(3.4)
dengan * menyatakan elemen diagonal yang didefisikan qii = −∑ qij j
b. Model Kimura Model ini diperkenalkan pada tahun 1980. Asumsi yang digunakan pada model sebelumnya, model Jukes dan Cantor, yaitu substitusi antara keempat basa nukleotid terjadi secara acak tidak selalu benar dalam praktiknya. Sebagai contoh, transisi (mutasi basa A menjadi basa G dan sebaliknya atau basa C
43
menjadi basa T dan sebaliknya) lebih banyak dari pada transversi (basa A menjadi basa C, basa A menjadi basa T, basa G menjadi basa C atau basa G menjadi basa T). Dalam model ini, digunakan dua parameter yaitu β untuk intensitas transisi dan γ untuk intensitas tranversi. Perhatikan diagram transisi berikut ini.
Dalam model ini, frekuensi basa pada saat setimbang diasumsikan sama
(π A = π C = π G = π T = 14 )
dan matriks intensitas substitusi didefinisikan
sebagai berikut ⎛* 1 κ 1⎞ ⎟ ⎜ ⎜1 * 1 κ⎟ Q=⎜ κ 1 * 1⎟ ⎟ ⎜ ⎜1 κ 1 *⎟ ⎠ ⎝ dengan κ =
(3.5)
β yaitu rasio antara intensitas transisi dan intensitas tranversi. γ
Jika intensitas tranversi lebih besar daripada intensitas transisi, maka κ → 0 . Dan sebaliknya jika intensitas transisi lebih besar daripada intensitas tranversi, maka κ → ∞ . c. Model Felsenstein Model ini adalah model yang diperkenalkan oleh Felsentein pada tahun 1981. Dalam model ini, frekuensi basa pada saat setimbang diasumsikan tidak sama dan tidak digunakan parameter intensitas substitusi. Matriks intensitas substitusi didefinisikan sebagai berikut,
44
⎛ * πC ⎜ * ⎜π Q=⎜ A π πC ⎜ A ⎜π ⎝ A πC
π G πT ⎞ ⎟ π G πT ⎟ * πT ⎟ ⎟ π G * ⎟⎠
(3.6)
dengan π A ≠ π C ≠ π G ≠ π T .
d. Model Hasegawa, Kishino, dan Yano Model ini diperkenalkan oleh Hasegawa, Kishino dan Yano pada tahun 1985. Model ini banyak digunakan karena dianggap lebih sesuai dengan kondisi sebenarnya (Siepel dan Haussler, 2003). Parameter yang digunakan adalah κ yaitu rasio antara intensitas transisi dan intensitas tranversi dan matriks intensitas substitusi didefinisikan sebagai berikut,
πC ⎛ * ⎜ * ⎜π Q=⎜ A κπ πC ⎜ A ⎜π ⎝ A κπ C
κπ G π T ⎞ ⎟ π G κπ T ⎟ * πT ⎟ ⎟ πG * ⎟⎠
(3.7)
dengan π A ≠ π C ≠ π G ≠ π T . Jadi secara umum, Table 1 memperlihatkan keempat model tersebut. Table 1 Model-model Filogenetik
45
Keterangan lebih jelas tentang model filogenetik dapat dilihat pada tesis Rifan (2007). Dalam tesis tersebut Rantai Markov kontinu digunakan untuk membangun pohon filogenetik khususnya untuk menghitung peluang perubahan basa yang terjadi dari spesies anchestral menjadi basa pada spesies saat ini. Ruang keadaan untuk Rantai Markov kontinu adalah basa nukleotid dalam barisan DNA yaitu basa A, C, G, dan T. Sedangkan indeks parameternya adalah waktu atau durasi perubahan basa dari spesies anchestral yang direpresentasikan dengan panjang cabang. Panjang cabang memiliki distribusi eksponensial sebagai akibat dari proses perubahan basa diasumsikan
mengikuti
proses
Poisson.
Pohon
filogenetik
dipilih
dengan
menggunakan metode Bayes. 3.5.
Likelihood Model Filogenetik
Seperti yang sudah dijelaskan sebelumnya, bahwa model filogenetik memiliki rate evolusi yang sama untuk setiap site. Hal ini mengakibatkan perhitungan panjang cabang pohon menjadi lepas dari rate evolusinya. Padahal rate evolusi ini bisa juga dikaitkan dengan panjang cabangnya. Untuk itu perlu diharapkan suatu pohon mempunyai rate evolusi yang berbeda di setiap site. Metode maksimum likelihood dapat digunakan untuk variasi rate evolusi. Untuk itu ada 2 hal yang harus diperhatikan, yaitu rate berbeda antar site dan terdapat korelasi antara rate evolusi pada site yang berdekatan. 46
Gambar 4 Representasi Model Filogenetik dan MMT
Dalam tesis ini akan dihitung likelihood filogenetik dengan rate evolusi yang mengikuti MMT (Gambar 4) dengan menggunakan algoritm Felsenstein. Rate evolusi diasumsikan berdistribusi gamma dengan parameter α dan β merupakan fungsi dari intensitas perubahan basa dan rasio antara intensitas transisi dengan tranversi. Sebenarnya asumsinya tidak harus gamma, bisa juga dengan menggunakan distribusi log-normal dan weibull. Ketiga distribusi di atas dipilih karena secara deskriptif mereka condong ke kiri (cenderung ke rate evolusi yang kecil). Hal ini cukup beralasan karena dalam filogenetik jarang sekali yang dibuat topologi pohonnya itu spesies-spesies yang memiliki kelas, ordo, atau famili yang berbeda. Jadi susunan barisan DNA-nya juga dominan yang sama. Selain itu, dalam susunan DNA, yang lebih dominan adalah intensitas transisi basanya dibandingkan dengan intensitas tranversinya. Untuk menentukan α dan β , model ini bisa dibagi ke dalam 2 kasus. Yang pertama dianggap tidak ada transisi, dan yang kedua tidak ada transisi dan transversi. Anggap rate dari 2 kasus ini adalah α dan β . Pada kasus I, karena tidak ada transisi, maka kita hanya cukup memperhatikan definisinya di purin dan pirimidin, berarti peluang basa A tetap menjadi A adalah peluangnya menjadi
πA , dan kalau diganti terhadap G, π A + πG
πG . Hal ini juga berlaku untuk basa C dan T. Pada kasus II, π A + πG
karena tidak ada transisi dan tranversi, berarti perubahan basanya memiliki 4
47
kemungkinan. Basa A diganti dengan A yang lain dengan peluang π A , diganti dengan G peluangnya menjadi π G , dan berlaku juga untuk C dan T. Dari kedua kasus tersebut, maka intensitas perubahan basa untuk setiap site dapat dituliskan sebagai berikut:
⎡
⎞ ⎛ π C ⎞⎤ ⎟ + πT ⎜ ⎟⎥ ⎠ ⎝ π C + π T ⎠⎦ ⎣ + β [π A (1 − π A ) + π G (1 − π G ) + π C (1 − π C ) + π T (1 − π T ) ] ⎛
⎛ πA ⎞ ⎛ πT πG ⎞ ⎟ + πG ⎜ ⎟ + πC ⎜ ⎝ π A + πG ⎠ ⎝ π A +πG ⎠ ⎝ πC + πT
θ = α ⎢π A ⎜
(3.8)
Misalkan π R dan π Y adalah proporsi perubahan basa dari purin dan pirimidin, sehingga kita bisa menuliskan
π R = π A +πG
πY = π C + πT
dan
(3.9)
maka persamaan (3.8) bisa disederhanakan menjadi
⎛ π Aπ G
θ = 2α ⎜
⎝ πR
+
π Cπ T ⎞ 2 2 2 2 ⎟ + β (1 − π A − π G − π C − π T ) πY ⎠
(3.10)
Rasio dari transisi dan tranversi bisa ditulis sebagai
⎛π π π π ⎞ 2α ⎜ A G + C T ⎟ + 2β (π Aπ G + π Cπ T ) πY ⎠ ⎝ πR κ= 2βπ Rπ Y
(3.11)
Dengan mensubtitusi persamaan (3.10) dan (3.11), maka diperoleh:
α=
2π Rπ Y κ − 2(π Aπ G + π Cπ T ) θ 1+ κ ⎛π π π π ⎞ 2⎜ A G + C T ⎟ πY ⎠ ⎝ πR
β=
1
(3.12)
θ
(3.13)
2π Rπ Y 1 + κ
Misalkan ci merupakan kategori rate yang terletak pada site i, dan ratenya sendiri didefinisikan sebagai rci , maka jika terdapat n site, kita bisa mendefinisikan likelihood dari model filogenetik ψ sebagai berikut: L = P (O ψ ) =
∑∑ ...∑ P ( c , c ,..., c ) P ( O ψ , r , r 1
c1
c2
cn
2
n
c1
c2
,..., rcn
)
(3.14)
48
n
=
∑∑ ...∑ P ( c , c ,..., c ) ∏ P ( O ψ , r ) 1
c1
c2
2
n
i
n
=
=
c2
i =1
cn
⎡
)
(3.16)
⎤ ∑ f ⎢∑ ...∑ P ( c ,..., c c ) ∏ P ( O ψ , r )⎥⎦
(3.17)
∑f
(3.18)
c1
c1
=
(
∑∑ ...∑ fc1 Pc1c2 Pc2c3 ...Pcn−1cn ∏ P Oi ψ , rci c1
(3.15)
ci
i =1
cn
c1
c1
n
⎣ c2
2
1
n
i
ci
i =1
cn
L(1) c1
Jika diambil sebanyak 3 kategori rate evolusi, dan terdapat 1000 site, maka akan ada sebanyak 31000 = 10477 kali penjumlahan. Hasil ini sangat tidak efisien dalam hal komputasi. Untuk itu, dengan menggunakan algoritma Felsenstein ini, operasi yang dibutuhkan akan jauh lebih mudah dan sedikit. Persamaan (3.14) berubah menjadi (3.5) dikarenakan adanya asumsi bahwa rate evolusi independen untuk setiap site, sehingga bisa dituliskan sebagai perkalian dari peluang masing-masing observasi bersyarat rate evolusi untuk site yang bersangkutan. Persamaan (3.15) menjadi (3.16) karena kategori rate berjalan dari sepanjang n site, sehingga bisa ditulis sebagai perkalian dari peluang prior c1 dengan matriks transisi antar kategori rate. Menuju ke persamaan (3.17) itu hanya merubah peluang gabungan menjadi peluang bersyarat, dan kemudian memisalkan […] menjadi L(1) c1 sehingga terbentuk persamaan (3.18). Sekarang kita definisikan L(1) c1 : L(1) = c1
n
(
∑ ...∑ P ( c2 ,..., cn c1 )∏ P Oi ψ , rci c2
i =1
cn
(
= P O1 ψ , rc1
n
) ∑ ...∑ P ( c ,..., c c )∏ P ( O ψ , r ) 2
c2
1
i
i =2
ci
n ⎡ ⎤ P ... P c ,..., c c ( 2 n 1 )∏ P Oi ψ , rci ⎥ ∑ c1c2 ⎢ ∑ ∑ c2 i =2 ⎣ c3 cn ⎦
(
)
(
)∑ P
= P O1 ψ , rc1
n
cn
(
= P O1 ψ , rc1
)
c2
c1c2
L(2) c2
)
(3.19)
Jika diperumum untiuk sebarang k, maka akan diperoleh
(
L(ckk ) = P Ok ψ , rck
)∑ P ck +1
ck ck +1
L(ckk ++11)
(3.20)
49
(
)
dan untuk k = n, L(cnn ) = P On ψ , rcn . Jadi langkah pertama yang dilakukan adalah
(
menentukan P Ok ψ , rck
)
untuk setiap k = 1, 2, ..., n yang tidak lain kalau di dalam
MMT lazim disebut dengan peluang emisi. Setelah itu menghitung L(ckk ) dari k = n, n1, ..., 1 sehingga terakhir kita dapat L(1) c1 . Terakhir, kita bisa menentukan likelihood filogenetik L dengan menggunakan persamaan (3.18). 3.6.
Kombinasi Rate Evolusi yang Optimal Setelah diperoleh likelihood dari filogenetik, satu hal yang mungkin menarik
perhatian kita adalah bagaimana variasi rate evolusi pada site yang berbeda. Likelihood telah dihitung dengan cara menjumlahkan semua peluang dari kombinasi yang mungkin dari rate evolusi. Dari semua kemungkinan tersebut, akan ada satu kombinasi yang memiliki kontribusi terbesar terhadap likelihood, dan kontribusi tersebut adalah:
(
R = max P ( c1 , c2 ,..., cn ) P O ψ , rc1 , rc2 ,..., rcn c1 ,c2 ,...,cn
)
(3.21)
Masalah ini bisa diselesaikan dengan menggunakan algoritma Viterbi. Analog dengan perhitungan L, kita juga bisa mendefinisikan Rc(kk ) sebagai berikut: Rc(kk )
{
(
= max P ( ck +1 ,..., cn ck ) P O ( k ) ψ , rck , rck +1 ,..., rcn ck +1 ,..., cn
(
)
{
= P Ok ψ , rck max Pck ck +1 Rc(kk++1 1) ck +1
}
)}
(3.22) (3.23)
Persamaan (3.23) analog dengan persamaan (3.20). Untuk k = n, definisi (3.13) menjadi
(
Rc(nn ) = P On ψ , rcn
)
(3.24)
Dengan memanfaatkan persamaan (3.23) dari n-1, n-2, dan terus ke bawah untuk semua kemungkinan kategori rate c1 pada sampai k = 1, akan diperoleh Rc(1) 1 site ke-1.
f c1 Rc(1) merupakan nilai yang memiliki kontribusi terbesar terhadap 1
likelihood dari semua kombinasi kategori rate. Meskipun R sudah bisa dihitung, tetapi kita masih belum mengetahui barisan dari kategori rate c1 , c2 ,..., cn . Akan tetapi, karena kita menggunakan persamaan (3.33) untuk setiap site, untuk setiap kemungkinan kategori rate di site tersebut yang
50
memiliki peluang maksimum, maka ck +1 akan menjadi kategori rate pada site berikutnya. Anggap peluang yang maksimum tersebut adalah Cc(kk ) , maka Cc(kk ) adalah nilai dari ck +1 yang terpilih melalui persamaan (3.23). Jika proses ini dilakukan terus sampai pada site ke-1, maka kita bisa mengetahui kategori rate c1.
Kita bisa
menggunakan Cc(1) untuk menentukan nilai c2 , dan kemudian Cc(2) digunakan untuk 1 2 mencari kategori rate pada site ke-3, dan begitu seterusnya. Setelah itu dilakukan backtracking sebagai bagian terakhir dari algoritma Viterbi. 3.7.
Estimasi Peluang Transisi Metode MMT digunakan karena adanya variasi rate evolusi untuk setiap site
yang berbeda. Salah satu parameter yang harus diestimasi adalah matriks transisi dari kategori rate. Dalam model ini, matriks transisi didefinisikan sebagai: Pij = λδ ij + (1 − λ ) f j
(3.25)
di mana λ adalah parameter otokorelasi. Setiap site diasumsikan memiliki peluang sebesar λ bahwa rate pada site tersebut adalah sama dengan rate pada site sebelumnya. Dan dengan peluang (1 − λ ) bahwa rate bersifat acak dari suatu distribusi yang setimbang. Dari metode ini, f j diindikasikan memiliki distribusi uniform. Pemilihan distribusi uniform ini karena dilihat dari cara kategori rate dipilih. Sedangkan δ ij merupakan fungsi delta Kronecker yang mana bernilai 1 untuk i = j dan 0 untuk yang lainnya. Jika λ bernilai mendekati 1, maka kategori rate di antara site yang bersebelahan memiliki otokorelasi yang besar, dan jika bernilai 0 maka tidak ada otokorelasi. Dalam studi kasus tidak akan dibahas tentang estimasi peluang transisi. Hasil dari model ini nantinya bisa digunakan untuk menentukan panjang cabang. Panjang cabang tersebut bisa dihitung dengan menggunakan metode NewtonRaphson (Felsenstein dan Churchill, 1996). Dan hasilnya bisa dibandingkan dengan metode standar seperti metode maksimum likelihood (Felsenstein, 1981). Namun hal ini tidak bisa kita bahas dalam tesis ini karena kedua metode di atas bisa cukup panjang untuk dikaji dan masalahnya akan menjadi lebih kompleks jika tetap dipaksakan untuk ditulis dalam tesis ini.
51
