Dosen Pengampu : Prof. Dr. Suhartono, Mkom

UIN Maulana Malik Ibrahim Malang

Mata Kuliah Bioinformatika

Dalam biologi molekuler, banyak situasi melibatkan peristiwa penghitungan: seperti berapa banyak kodon (kode genetik) yang menggunakan ejaan tertentu, berapa banyak pembacaan DNA yang cocok dengan referensi, berapa banyak diagram CG yang diamati dalam urutan DNA. Hitungan ini memberi kita variabel diskrit , berbeda dengan kuantitas seperti massa dan intensitas yang diukur pada skala kontinu.

Jika kita ingin mengetahui seberapa sering 3 mutasi dapat terjadi pada model Poisson (5), kita dapat menggunakan fungsi R untuk menghasilkan probabilitas. Dimana x=3 peristiwa dan nilai parameter laju distribusi Poisson, yang disebut lambda menjadi 5 .

dpois(x = 3, lambda = 5)
## [1] 0.1403739

Artinya peluang untuk melihat tepat 3 peristiwa adalah sekitar 0,14

Jika kita ingin membangkitkan probabilitas semua nilai dari 0 hingga 12, kita tidak perlu menulis perulangan. Kita cukup mengatur argumen pertama menjadi vektor dari 13 nilai ini, menggunakan operator urutan R, titik dua “ :”. Kita dapat melihat probabilitas dengan memplot datanya.

0:12
##  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12
dpois(x = 0:12, lambda = 5)
##  [1] 0.006737947 0.033689735 0.084224337 0.140373896 0.175467370 0.175467370
##  [7] 0.146222808 0.104444863 0.065278039 0.036265577 0.018132789 0.008242177
## [13] 0.003434240
barplot(dpois(0:12, 5), names.arg = 0:12, col = "blue")

Grafik diatas menunjukan probabilitas munculnya jumlah mutasi 0 hingga 12, seperti yang dimodelkan oleh distribusi Poisson (5). Berdasarkan plot histogram diatas,menunjukan bahwa mutasi sebanyak 4 atau 5 akan sering muncul atau terlihat dan mutasi sebanyak 12 atau lebih akan jarang muncul atau terlihat.Distribusi Poisson adalah model yang baik untuk kejadian langka seperti mutasi.

Menggunakan Model Probabilitas Diskrit

Tidak semua peristiwa adalah biner. Misalnya, genotipe dalam organisme diploid dapat mengambil tiga tingkatan (AA, Aa, aa).Ketika kita mengukur variabel kategori pada sampel, kita sering ingin menghitung frekuensi tingkat yang berbeda dalam vektor hitungan. R memiliki pengkodean khusus untuk variabel kategori dan menyebutnya faktor 6 . Di sini kami menangkap genotipe darah yang berbeda untuk 19 subjek dalam vektor yang kami tabulasi

genotype = c("AA","AO","BB","AO","OO","AO","AA","BO","BO",
             "AO","BB","AO","BO","AB","OO","AB","BB","AO","AO")
table(genotype)
## genotype
## AA AB AO BB BO OO 
##  2  2  7  3  3  2

Saat membuat faktor , R secara otomatis mendeteksi level. Anda dapat mengakses level dengan levels fungsi.

genotypeF = factor(genotype)
levels(genotypeF)
## [1] "AA" "AB" "AO" "BB" "BO" "OO"
table(genotypeF)
## genotypeF
## AA AB AO BB BO OO 
##  2  2  7  3  3  2

1. Percobaan Bernoulli

Melempar koin memiliki dua kemungkinan hasil. Eksperimen sederhana ini, yang disebut percobaan Bernoulli, dimodelkan menggunakan variabel acak Bernoulli.

Beberapa eksperimen untuk melihat variabel acak menggunakan fungsi R khusus yang disesuaikan untuk menghasilkan setiap jenis distribusi. Semuanya dimulai dengan huruf r, diikuti dengan spesifikasi model, di sini menggunakan rbinom, di mana binom singkatan yang digunakan untuk binomial.

Contoh mensimulasikan urutan 15 lemparan koin yang sama. Untuk mendapatkan hasil dari 15 percobaan Bernoulli dengan probabilitas keberhasilan sama dengan 0,5 (koin yang sama).

rbinom(15, prob = 0.5, size = 1)
##  [1] 0 1 1 1 1 1 0 0 1 1 1 0 1 0 1

Contoh mensimulasikan 12 percobaan melempar bola ke dalam 2 kotak seperti yang ditunjukkan pada Gambar dibawah ini , dengan probabilitas jatuh di kotak sebelah kanan 2/3 dan di kotak sebelah kiri 1/3.

knitr::include_graphics("gambar_bernauli.png")

rbinom(12, prob = 2/3, size = 1)
##  [1] 1 1 1 1 0 1 1 1 1 0 0 1

PERHITUNGAN KEBERHASILAN BINOMIAL

Jika yang diperhatikan hanya jumlah bola yang masuk ke kotak sebelah kanan, maka urutan lemparan tidak masalah dan kita bisa mendapatkan angka ini hanya dengan mengambil jumlah sel dalam vektor keluaran.

rbinom(1, prob = 2/3, size = 12)
## [1] 7

Banyaknya keberhasilan dalam 15 percobaan Bernoulli dengan probabilitas keberhasilan 0,3 disebut variabel acak binomial atau variabel acak yang mengikuti distribusi.

set.seed(235569515)
rbinom(1, prob = 0.3, size = 15)
## [1] 5

Distribusi massa probabilitas lengkap tersedia dengan mengetik:

probabilities = dbinom(0:15, prob = 0.3, size = 15)
round(probabilities, 2)
##  [1] 0.00 0.03 0.09 0.17 0.22 0.21 0.15 0.08 0.03 0.01 0.00 0.00 0.00 0.00 0.00
## [16] 0.00

Tampilann grafik Histogram dari hasil distribusi diatas :

barplot(probabilities, names.arg = 0:15, col = "green")

2. Distribusi Poisson

Ketika peluang sukses kecil dan jumlah percobaan n besar, distribusi binomial B (n,p) dapat didekati dengan tepat dengan distribusi yang lebih sederhana, yaitu distribusi Poisson dengan parameter laju (lamda = np).

Contoh Simulasi proses mutasi sepanjang 10.000 posisi dengan tingkat mutasi (5 x 10^-4) dan menghitung jumlah mutasinya. Kemudian diplot distribusinya dengan bar plot fungsi.

rbinom(1, prob = 5e-4, size = 10000)
## [1] 6
simulations = rbinom(n = 300000, prob = 5e-4, size = 10000)
barplot(table(simulations), col = "red")

MODEL GENERATIF UNTUK DETEKSI EPITOP

Epitop merupakan bagian spesifik dari antigen makromolekul yang mengikat antibodi. Dalam kasus antigen protein yang dikenali oleh sel T, epitop atau determinannya adalah bagian atau situs peptida yang berikatan dengan molekul Major Histocompatibility Complex (MHC) untuk dikenali oleh reseptor sel T (TCR).

knitr::include_graphics("gambar_antibodi.png")

Antibodi adalah sejenis protein yang dibuat oleh sel darah putih tertentu sebagai respons terhadap zat asing di dalam tubuh, yang disebut antigen.Antibodi mengikat (dengan lebih atau kurang spesifisitas) ke antigennya. Tujuan dari pengikatan adalah untuk membantu menghancurkan antigen. Antibodi dapat bekerja dalam beberapa cara, tergantung pada sifat antigen. Beberapa antibodi menghancurkan antigen secara langsung.

Contoh simulasi 50 variabel bebas Bernoulli dengan p = 0,01 dengan pendekatan distribusi Poisson (parameter variabel acak 0,5).

load("e100.RData")
barplot(e100, ylim = c(0, 7), width = 0.7, xlim = c(-0.5, 100.5),
  names.arg = seq(along = e100), col = "darkolivegreen")

Untuk melihat peluang angka 7 (atau lebih besar) dengan mempertimbangkan Poisson (0,5) variabel acak dapat dihitung dalam bentuk tertutup, fungsi distribusi kumulatif.

1 - ppois(6, 0.5)
## [1] 1.00238e-06
ppois(6, 0.5, lower.tail = FALSE)
## [1] 1.00238e-06

MENGHITUNG PROBABILITAS DENGAN SIMULASI

Perhitungan probabilitas teoritis cukup sederhana dan dapat diketahui hasilnya dengan perhitungan eksplisit. Dalam praktiknya, kadang cenderung lebih rumit, dan lebih baik menghitung probabilitas menggunakan metode Monte Carlo, yaitu simulasi komputer berdasarkan model generatif,dimana kita yang menemukan probabilitas kejadian yang kita minati. Berikut adalah contohnya, dimana kita menghasilkan 100.000 contoh dan memilih maksimum dari 100 nomor terdistribusi Poisson.

maxes = replicate(100000, {max(rpois(100, 0.5))})
table(maxes)
## maxes
##     1     2     3     4     5     6     7     9 
##     7 23027 60840 14365  1604   141    15     1

Jika dalam 16 dari 100000 percobaan, Probabilitas nilai maksimum lebih dari sama dengan 7 adalah

mean( maxes >= 7 )
## [1] 0.00016

3. Distribusi multinomial: kasus DNA

Jika suatu model mempunyai empat kemungkinan hasil, seperti kotak pada Gambar dibawah ini, atau ketika mempelajari jumlah empat nukleotida A,C,G dan T , kita perlu memperluas model binomial.Sebagai contoh ketika melempar bola ke dalam kotak dengan ukuran berbeda yang sesuai dengan probabilitas yang berbeda, kita dapat memberi label probabilitas ini pA, pB, pG, pT. Sama seperti dalam kasus binomial jumlah probabilitas semua hasil yang mungkin adalah 1, pA + pB + pG + pT = 1

knitr::include_graphics("gambar_multinomial.png")

Contoh soal :

Misalkan kita memiliki empat kotak yang kemungkinannya sama. Dengan menggunakan rumus, berapakah peluang untuk melihat 4 di kotak pertama, 2 di kotak kedua, dan tidak ada di dua kotak lainnya?

Jawab :

Kita sering menjalankan eksperimen simulasi untuk memeriksa apakah data yang kita lihat konsisten dengan kemungkinan model empat kotak paling sederhana di mana setiap kotak memiliki probabilitas 1/4 yang sama. Dalam beberapa hal itu adalah strawman (tidak ada yang menarik yang terjadi). Kita akan melihat lebih banyak contoh tentang ini di Bab 2 . Di sini kita menggunakan beberapa perintah R untuk menghasilkan vektor hitungan seperti itu. Pertama misalkan kita memiliki 8 karakter dari empat jenis yang berbeda, kemungkinan yang sama:

pvec = rep(1/4, 4)
t(rmultinom(1, prob = pvec, size = 8))
##      [,1] [,2] [,3] [,4]
## [1,]    2    2    3    1
rmultinom(n = 8, prob = pvec, size = 1)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]    1    0    1    0    0    0    0    0
## [2,]    0    0    0    1    0    0    1    0
## [3,]    0    1    0    0    1    0    0    1
## [4,]    0    0    0    0    0    1    0    0
rmultinom(n = 1, prob = pvec, size = 8)
##      [,1]
## [1,]    2
## [2,]    2
## [3,]    4
## [4,]    0

Contoh dalam 1000 simulasi dari hipotesis nol menggunakan fungsi rmultinomial, data ditampilkan dalam 11 kolom pertama. Setiap kolom dalam matriks obsunder0 adalah contoh simulasi. Angka-angka dalam kotak sangat variasi, nilai yang diharapkan adalah 5 = 20/4.

obsunder0 = rmultinom(1000, prob = pvec, size = 20)
dim(obsunder0)
## [1]    4 1000
obsunder0[, 1:11]
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,]    4    6    5    4    4    6    3    6    4     5     3
## [2,]    3    2    6    4    5    6    4    3    8     4     5
## [3,]    8    7    6    5    6    5    6    6    3     5     5
## [4,]    5    5    3    7    5    3    7    5    5     6     7
expected0 = pvec * 20
sum((obsunder0[, 1] - expected0)^2 / expected0)
## [1] 2.8
sum((obsunder0[, 2] - expected0)^2 / expected0)
## [1] 2.8
sum((obsunder0[, 3] - expected0)^2 / expected0)
## [1] 1.2
stat = function(obsvd, exptd = 20 * pvec) {
  sum((obsvd - exptd)^2 / exptd)
}
stat(obsunder0[, 1])
## [1] 2.8
S0 = apply(obsunder0, 2, stat)
summary(S0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.200   2.800   3.086   4.400  18.800
hist(S0, breaks = 25, col = "lavender", main = "")

q95 = quantile(S0, probs = 0.95)
q95
## 95% 
## 7.6