Dosen Pengampu : Prof. Dr. Suhartono, Mkom

UIN Maulana Malik Ibrahim Malang

Mata Kuliah Bioinformatika

Dalam banyak situasi nyata, baik model generatif maupun parameter yang tidak diketahui, kita perlu memperkirakan dengan menggunakan data yang telah dikumpulkan. Pemodelan statistik bekerja dari data ke atas (upwards) ke model yang mungkin menjelaskan data secara masuk akal (plausibly). Langkah penalaran ke atas ini disebut inferensi statistik . Bab ini akan menunjukkan kepada kita beberapa distribusi dan mekanisme estimasi yang berfungsi sebagai blok bangunan untuk inferensi. Meskipun contoh dalam bab ini semuanya parametrik (yaitu, model statistik hanya memiliki sejumlah kecil parameter yang tidak diketahui), prinsip-prinsip yang kita diskusikan akan digeneralisasi.

Contoh sederhana pemodelan statistik

load("e100.RData")
e99 = e100[-which.max(e100)]
e99
##  [1] 2 0 1 0 0 0 2 1 0 0 0 0 0 0 0 0 0 1 1 1 0 0 1 0 1 0 0 0 1 0 0 0 0 1 0 1 1 0
## [39] 1 2 2 1 0 2 0 1 0 1 1 1 1 0 1 0 0 0 0 1 2 2 1 0 0 0 0 0 1 0 0 0 1 1 0 0 1 1
## [77] 0 0 0 0 1 0 0 0 1 0 1 0 0 1 1 0 0 1 1 0 0 1 0

Distribusi data epitop yang diamati tanpa outlier

barplot(table(e99), space = 0.8, col = "chartreuse4")

Satu diagram kesesuaian visual dikenal sebagai rootogram ( William S. Cleveland 1988 ) ; Diagram batang yang menggantung dengan jumlah yang diamati dari titik merah teoretis. Jika hitungan sesuai persis dengan nilai teoretisnya, bagian bawah kotak akan sejajar persis dengan sumbu horizontal.

library("vcd")
## Loading required package: grid
gf1 = goodfit( e99, "poisson")
rootogram(gf1, xlab = "", rect_gp = gpar(fill = "chartreuse4"))

Mengestimasi parameter distribusi Poisson

Berapa nilai rata-rata Poisson yang membuat data paling mungkin terjadi? Pada langkah pertama, menghitung hasilnya.

table(e100)
## e100
##  0  1  2  7 
## 58 34  7  1

Kemudian mencoba nilai yang berbeda untuk mean Poisson dan melihat mana yang paling cocok dengan data kita. Jika maksudnya lamda dari distribusi Poisson adalah 3, hitungannya akan terlihat seperti ini :

databaru <- rpois(100, 3)
databaru
##   [1] 2 6 6 3 1 1 1 2 2 0 3 1 1 2 1 2 4 3 4 4 6 1 3 0 6 4 1 6 5 3 0 2 1 5 1 2 4
##  [38] 4 3 3 4 5 4 3 4 3 1 1 0 3 3 1 4 3 6 4 4 1 2 2 1 3 3 4 5 2 0 0 3 5 5 0 1 4
##  [75] 3 1 4 3 8 4 3 4 5 5 2 5 4 3 2 2 1 7 1 4 3 1 5 0 4 5

Dari tabel diatas, terlihat data 2 dan 3 lebih sering muncul. JIka kita melihat lamda =3 tidak mungkin menghasilkan data ini, karena penghitungannya tidak begitu cocok.

table(databaru)
## databaru
##  0  1  2  3  4  5  6  7  8 
##  8 20 13 20 20 11  6  1  1

Menghitung probabilitas melihat data jika nilai parameter Poisson adalah m. Karena kami menganggap data berasal dari penarikan independen, probabilitas ini hanyalah produk dari probabilitas individu

prod(dpois(c(0, 1, 2, 7), lambda = 3) ^ (c(58, 34, 7, 1)))
## [1] 1.392143e-110

Menghitung kemungkinan untuk banyak nilai parameter Poisson yang berbeda. Untuk melakukan ini, kita perlu menulis fungsi kecil yang menghitung probabilitas data untuk nilai yang berbeda.

loglikelihood  =  function(lambda, data = e100) {
  sum(log(dpois(data, lambda)))
}

Sekarang kita dapat menghitung kemungkinan untuk seluruh rangkaian lambdanilai dari 0,05 hingga 0,95

lambdas = seq(0.05, 0.95, length = 100)
loglik = vapply(lambdas, loglikelihood, numeric(1))
plot(lambdas, loglik, type = "l", col = "red", ylab = "", lwd = 2,
     xlab = expression(lambda))
m0 = mean(e100)
abline(v = m0, col = "blue", lwd = 2)
abline(h = loglikelihood(m0), col = "purple", lwd = 2)

m0
## [1] 0.55

Fungsi vapply yang dilakukan dalam kode di atas :

Menggunakan fungsi goodfit sebagai jalan singkatnya.

gf  =  goodfit(e100, "poisson")
names(gf) 
## [1] "observed" "count"    "fitted"   "type"     "method"   "df"       "par"
gf$par
## $lambda
## [1] 0.55

Output dari goodfitadalah objek komposit yang disebut daftar. Salah satu komponennya disebut pardan berisi nilai parameter yang dipasang untuk distribusi yang dipelajari. Dalam hal ini hanya satu angka, perkiraan.

gf$observed
## [1] 58 34  7  0  0  0  0  1
gf$fitted
## [1] 5.769498e+01 3.173224e+01 8.726366e+00 1.599834e+00 2.199771e-01
## [6] 2.419749e-02 2.218103e-03 1.742795e-04
gf$type
## [1] "poisson"
cb <- rbinom(120, prob = 0.8, size = 1)
cb
##   [1] 0 0 1 1 0 0 1 1 1 1 0 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1
##  [38] 0 1 0 1 1 1 0 1 1 0 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
##  [75] 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 0 0 0 1 1 1 1 1 1 1 0 0 1 1 1 1 0
## [112] 1 1 1 1 0 1 0 0 0

Distribusi binomial dan kemungkinan maksimum

Dalam distribusi binomial ada dua parameter: jumlah percobaan, yang biasanya diketahui, dan probabilitasmelihat 1 dalam percobaan. Probabilitas ini sering tidak diketahui.

Sebagai contoh, kita mengambil sampel darijantan dan menguji mereka untuk buta warna merah-hijau. Kita dapat mengkodekan data sebagai 0 jika subjek tidak buta warna dan 1 jika dia buta warna. Kami meringkas data dengan tabel :

table(cb)
## cb
##  0  1 
## 29 91

Dalam kasus khusus ini, intuisi kita mungkin akan memberi perkiraan, yang ternyata merupakan perkiraan kemungkinan maksimum. Kami meletakkan topi di atas surat itu untuk mengingatkan kami bahwa ini bukan (harus) nilai sebenarnya yang mendasarinya, tetapi perkiraan yang kami buat dari data.

Seperti sebelumnya dalam kasus Poisson, jika kita menghitung kemungkinan untuk banyak kemungkinan, kita dapat memplotnya dan melihat di mana jatuh maksimumnya

probs  =  seq(0, 0.3, by = 0.005)
likelihood = dbinom(sum(cb), prob = probs, size = length(cb))
plot(probs, likelihood, pch = 16, xlab = "probability of success",
       ylab = "likelihood", cex=0.6)

probs[which.max(likelihood)]
## [1] 0.3

Kemungkinan untuk distribusi binomial

emungkinan dan probabilitas adalah fungsi matematis yang sama, hanya ditafsirkan dengan cara yang berbeda – dalam satu kasus, fungsi tersebut memberi tahu kita seberapa besar kemungkinan untuk melihat sekumpulan nilai data tertentu, mengingat parameternya; dalam kasus lain, kami menganggap data seperti yang diberikan, dan meminta nilai parameter yang kemungkinan menghasilkan data ini. Memperkirakan n = 300 dan kita amati y = 40 keberhasilan. Kemudian, untuk distribusi binomial menggunakan fungsi:

loglikelihood = function(p, n = 300, y = 40) {
  log(choose(n, y)) + y * log(p) + (n - y) * log(1 - p)
}

Yang diplot untuk kisaran probabilitas (p) dari 0 hingga 1

p_seq = seq(0, 1, by = 0.001)
plot(p_seq, loglikelihood(p_seq), xlab = "p", ylab = "log f(p|y)", type = "l")

Pemodelan jumlah DNA: pasangan basa

Ada empat molekul dasar DNA: A - adenin, C - sitosin, G - guanin, T - timin. Nukleotida diklasifikasikan menjadi 2 kelompok: purin ( A dan G ) dan pirimidin ( C dan T ). Binomial akan berfungsi sebagai model untuk pengelompokan purin/pirimidin tetapi tidak jika kita ingin menggunakan A, C, G, T ; untuk itu kita membutuhkan model multinomial dari 1.4 . Mari kita lihat pola nyata yang terjadi pada frekuensi ini.

Bias nukleotida

Bagian ini menggabungkan estimasi dan pengujian dengan simulasi dalam contoh nyata. Data dari satu untai DNA untuk gen bakteri Staphylococcus aureus tersedia dalam file fastastaphsequence.ffn.txt , yang dapat kita baca dengan fungsi dari paket Biokonduktor Biostrings.

library("Biostrings")
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
## Loading required package: stats4
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## The following object is masked from 'package:vcd':
## 
##     tile
## The following object is masked from 'package:grDevices':
## 
##     windows
## Loading required package: XVector
## Loading required package: GenomeInfoDb
## 
## Attaching package: 'Biostrings'
## The following object is masked from 'package:grid':
## 
##     pattern
## The following object is masked from 'package:base':
## 
##     strsplit
staph = readDNAStringSet("staphsequence.ffn.txt", "fasta")

Melihat gen pertama :

staph[1]
## DNAStringSet object of length 1:
##     width seq                                               names               
## [1]  1362 ATGTCGGAAAAAGAAATTTGGGA...AAAAAGAAATAAGAAATGTATAA lcl|NC_002952.2_c...
length(staph)
## [1] 2650
letterFrequency(staph[[1]], letters = "ACGT", OR = 0)
##   A   C   G   T 
## 522 219 229 392

Menguji apakah nukleotida terdistribusi secara merata di empat nukleotida untuk gen pertama. Karena sifat fisiknya yang berbeda, seleksi evolusioner dapat bekerja pada frekuensi nukleotida. Jadi kita bisa bertanya apakah, katakanlah, sepuluh gen pertama dari data ini berasal dari multinomial yang sama. Kami tidak memiliki referensi sebelumnya, kami hanya ingin memutuskan apakah nukleotida terjadi dalam proporsi yang sama pada 10 gen pertama. Jika tidak, ini akan memberi kita bukti untuk berbagai tekanan selektif pada sepuluh gen ini.

letterFrq = vapply(staph, letterFrequency, FUN.VALUE = numeric(4),
         letters = "ACGT", OR = 0)
colnames(letterFrq) = paste0("gene", seq(along = staph))
tab10 = letterFrq[, 1:10]
computeProportions = function(x) { x/sum(x) }
prop10 = apply(tab10, 2, computeProportions)
round(prop10, digits = 2)
##   gene1 gene2 gene3 gene4 gene5 gene6 gene7 gene8 gene9 gene10
## A  0.38  0.36  0.35  0.37  0.35  0.33  0.33  0.34  0.38   0.27
## C  0.16  0.16  0.13  0.15  0.15  0.15  0.16  0.16  0.14   0.16
## G  0.17  0.17  0.23  0.19  0.22  0.22  0.20  0.21  0.20   0.20
## T  0.29  0.31  0.30  0.29  0.27  0.30  0.30  0.29  0.28   0.36
p0 = rowMeans(prop10)
p0
##         A         C         G         T 
## 0.3470531 0.1518313 0.2011442 0.2999714

misalkan p0 adalah vektor probabilitas multinomial untuk semua sepuluh gen dan menggunakan simulasi Monte Carlo untuk menguji apakah penyimpangan antara frekuensi huruf yang diamati dan nilai yang diharapkan di bawah anggapan ini berada dalam kisaran yang masuk akal.

Kami menghitung jumlah yang diharapkan dengan mengambil outer produk dari vektor probabilitas p0 dengan jumlah jumlah nukleotida dari masing-masing 10 kolom.

cs = colSums(tab10)
cs
##  gene1  gene2  gene3  gene4  gene5  gene6  gene7  gene8  gene9 gene10 
##   1362   1134    246   1113   1932   2661    831   1515   1287    696
expectedtab10 = outer(p0, cs, FUN = "*")
round(expectedtab10)
##   gene1 gene2 gene3 gene4 gene5 gene6 gene7 gene8 gene9 gene10
## A   473   394    85   386   671   924   288   526   447    242
## C   207   172    37   169   293   404   126   230   195    106
## G   274   228    49   224   389   535   167   305   259    140
## T   409   340    74   334   580   798   249   454   386    209

Membuat tabel acak dengan jumlah kolom yang benar menggunakan rmultinomfungsi. Tabel ini dihasilkan sesuai dengan hipotesis nol bahwa proporsi yang benar diberikan oleh p0.

randomtab10 = sapply(cs, function(s) { rmultinom(1, s, p0) } )
all(colSums(randomtab10) == cs)
## [1] TRUE

Sekarang kita ulangi B = 1000 kali. Untuk setiap tabel, menghitung statistik pengujian (fungsi stat) dan menyimpan hasilnya dalam vektor simulstat. Bersama-sama, nilai-nilai ini membentuk distribusi nol kami, karena mereka dihasilkan di bawah hipotesis nol yang p0merupakan vektor proporsi multinomial untuk masing-masing dari 10 gen.

stat = function(obsvd, exptd = 20 * pvec) {
   sum((obsvd - exptd)^2 / exptd)
}
B = 1000
simulstat = replicate(B, {
  randomtab10 = sapply(cs, function(s) { rmultinom(1, s, p0) })
  stat(randomtab10, expectedtab10)
})
S1 = stat(tab10, expectedtab10)
sum(simulstat >= S1)
## [1] 0
hist(simulstat, col = "lavender", breaks = seq(0, 75, length.out=50))
abline(v = S1, col = "red")
abline(v = quantile(simulstat, probs = c(0.95, 0.99)),
       col = c("darkgreen", "blue"), lty = 2)

Gambar diatas merupakan Histogram dari simulstat. Nilai dari S1ditandai dengan garis merah vertikal, nilai dari kuantil 0,95 dan 0,99 (lihat bagian berikutnya) dengan garis putus-putus.Kemungkinan melihat nilai sebesar S1=70,1 sangat kecil di bawah model nol . Itu terjadi 0 kali dalam 1000 simulasi kami bahwa nilai sebesar S1terjadi. Jadi kesepuluh gen tersebut tampaknya tidak berasal dari model multinomial yang sama.

Kuantil dan plot kuantil

Membandingkan simulstat nilai dan 1000 yang dihasilkan secara acak, dengan menampilkannya dalam histogram dengan masing-masing 50 bin. Kemudian menghitung nilai kuantil dan nilai simulsta.

qs = ppoints(100)
quantile(simulstat, qs)
##     0.5%     1.5%     2.5%     3.5%     4.5%     5.5%     6.5%     7.5% 
## 14.55860 15.54484 16.84526 17.80901 18.16761 18.61371 19.06581 19.60380 
##     8.5%     9.5%    10.5%    11.5%    12.5%    13.5%    14.5%    15.5% 
## 19.92724 20.34075 20.72252 20.88896 21.32175 21.57227 21.89892 22.09682 
##    16.5%    17.5%    18.5%    19.5%    20.5%    21.5%    22.5%    23.5% 
## 22.35680 22.57258 22.73162 22.91569 23.11727 23.43601 23.65418 23.85763 
##    24.5%    25.5%    26.5%    27.5%    28.5%    29.5%    30.5%    31.5% 
## 24.09073 24.27657 24.62778 24.76164 24.86827 25.15513 25.29018 25.47002 
##    32.5%    33.5%    34.5%    35.5%    36.5%    37.5%    38.5%    39.5% 
## 25.77724 25.98303 26.08868 26.38377 26.55398 26.64685 26.89372 27.08754 
##    40.5%    41.5%    42.5%    43.5%    44.5%    45.5%    46.5%    47.5% 
## 27.35285 27.59963 27.81148 27.96591 28.14593 28.42354 28.63637 28.82854 
##    48.5%    49.5%    50.5%    51.5%    52.5%    53.5%    54.5%    55.5% 
## 29.02865 29.12638 29.32464 29.56604 29.79535 30.03829 30.33925 30.65696 
##    56.5%    57.5%    58.5%    59.5%    60.5%    61.5%    62.5%    63.5% 
## 30.80166 30.96427 31.11338 31.27577 31.46880 31.78995 32.12662 32.43476 
##    64.5%    65.5%    66.5%    67.5%    68.5%    69.5%    70.5%    71.5% 
## 32.62955 32.77037 33.00507 33.23821 33.48136 33.90551 34.19280 34.48902 
##    72.5%    73.5%    74.5%    75.5%    76.5%    77.5%    78.5%    79.5% 
## 34.76996 35.24804 35.49960 35.73564 35.95940 36.08208 36.25534 36.55159 
##    80.5%    81.5%    82.5%    83.5%    84.5%    85.5%    86.5%    87.5% 
## 36.77663 37.02084 37.43979 37.89282 38.14887 38.56059 38.87871 39.18855 
##    88.5%    89.5%    90.5%    91.5%    92.5%    93.5%    94.5%    95.5% 
## 39.78180 40.35171 41.04222 41.63916 42.68649 43.07824 43.96040 45.04175 
##    96.5%    97.5%    98.5%    99.5% 
## 46.00812 46.62807 48.34254 53.25526
quantile(qchisq(qs, df = 30), qs)
##     0.5%     1.5%     2.5%     3.5%     4.5%     5.5%     6.5%     7.5% 
## 14.74308 16.23869 17.16382 17.87179 18.45879 18.96730 19.42030 19.83175 
##     8.5%     9.5%    10.5%    11.5%    12.5%    13.5%    14.5%    15.5% 
## 20.21086 20.56401 20.89588 21.20996 21.50896 21.79501 22.06984 22.33486 
##    16.5%    17.5%    18.5%    19.5%    20.5%    21.5%    22.5%    23.5% 
## 22.59125 22.83999 23.08192 23.31776 23.54813 23.77359 23.99461 24.21163 
##    24.5%    25.5%    26.5%    27.5%    28.5%    29.5%    30.5%    31.5% 
## 24.42502 24.63512 24.84224 25.04668 25.24867 25.44846 25.64627 25.84230 
##    32.5%    33.5%    34.5%    35.5%    36.5%    37.5%    38.5%    39.5% 
## 26.03673 26.22975 26.42152 26.61219 26.80192 26.99084 27.17910 27.36683 
##    40.5%    41.5%    42.5%    43.5%    44.5%    45.5%    46.5%    47.5% 
## 27.55414 27.74118 27.92804 28.11486 28.30175 28.48883 28.67621 28.86400 
##    48.5%    49.5%    50.5%    51.5%    52.5%    53.5%    54.5%    55.5% 
## 29.05231 29.24127 29.43100 29.62161 29.81322 30.00595 30.19994 30.39530 
##    56.5%    57.5%    58.5%    59.5%    60.5%    61.5%    62.5%    63.5% 
## 30.59219 30.79073 30.99107 31.19337 31.39779 31.60450 31.81367 32.02550 
##    64.5%    65.5%    66.5%    67.5%    68.5%    69.5%    70.5%    71.5% 
## 32.24019 32.45796 32.67904 32.90367 33.13213 33.36471 33.60172 33.84351 
##    72.5%    73.5%    74.5%    75.5%    76.5%    77.5%    78.5%    79.5% 
## 34.09046 34.34299 34.60155 34.86665 35.13886 35.41883 35.70725 36.00496 
##    80.5%    81.5%    82.5%    83.5%    84.5%    85.5%    86.5%    87.5% 
## 36.31286 36.63203 36.96368 37.30925 37.67041 38.04916 38.44789 38.86951 
##    88.5%    89.5%    90.5%    91.5%    92.5%    93.5%    94.5%    95.5% 
## 39.31761 39.79669 40.31253 40.87266 41.48728 42.17057 42.94329 43.83752 
##    96.5%    97.5%    98.5%    99.5% 
## 44.90723 46.25504 48.12235 51.45778

Memplot kuantil simulstat nilai, yang kami simulasikan di bawah hipotesis nol, terhadap distribusi nol teoretis.

qqplot(qchisq(ppoints(B), df = 30), simulstat, main = "",
  xlab = expression(chi[nu==30]^2), asp = 1, cex = 0.5, pch = 16)
abline(a = 0, b = 1, col = "red")