Gibbs Sampling

Pada post ini, saya akan menjelaskan mengenai implementasi algoritma Gibbs sampling untuk mendeteksi pola pada deret DNA atau populer dengan istilah motif finding seperti yang dijabarkan oleh Lawrence di paper-nya pada tahun 1993 Detecting Subtle Sequence Signals: A Gibbs Sampling Strategy for Multiple Alignment.

Sekilas mengenai Gibbs sampling. Gibbs sampling adalah salah satu algoritma keluarga dari Markov Chain Monte Carlo (MCMC). Kurang lebih intinya adalah kita bisa menghitung joint probability distribution dengan cara melakukan sampling satu per satu terhadap setiap variabel dengan berdasarkan nilai variabel lainnya alias full conditional probability. Materi pengenalan Gibbs sampling yang lebih lanjut bisa dipelajari di sini dan di sini.

Motif itu sendiri adalah pola yang terdapat pada beberapa deret DNA yang memiliki tingkat kemiripan yang tinggi. Contohnya bisa dilihat pada gambar di bawah. Seperti yang bisa dilihat, terdapat kemunculan pola dengan bentuk CAGGAT pada 3 deret dengan lokasi yang berbeda-beda. Motif finding bertujuan menemukan pola tersebut berdasarkan data semua deret string yang dimiliki. Karakter yang muncul selain pada motif biasa disebut background.

motif

Nah, bagaimana caranya kita bisa menggunakan Gibbs sampling untuk melakukan deteksi motif? Ada 2 asumsi sebelum kita memulai memakan Gibbs sampling. Pertama, pada satu deret DNA hanya ada 1 buah pola dan kedua, kita mengetahui berapa panjang dari pola tersebut. Kita akan implementasi Gibbs sampling dengan menggunakan bahasa pemrograman Python.

Pertama, kita inisialisasi state posisi secara acak

pos = [randint(0, len(seq)-w+1) for seq in sequences]

Kemudian kita lakukan iterasi sampai posisi yang kita cari konvergen. Kita lakukan iterasi terhadap semua deret satu per satu. Sebutlah deret yang aktif pada iterasi ini sebagai i.

    last_pos = None
    while last_pos != pos:
        last_pos = pos[:]
        
        # We pick the sequence, well, in sequence starting from index 0
        for i in xrange(K):

Kemudian dari semua deret yang kita punya, selain deret yang sedang aktif dalam iterasi, kita hitung model probabilitas motif dan background. Model untuk background cukup dihitung dengan menghitung berapa kali karakter suatu karakter muncul. Sedangkan untuk model motif, kita harus menghitung berapa kali sebuah karakter muncul pada posisi ke j pada motif tersebut.

 seq_minus = sequences[:]; del seq_minus[i]
 pos_minus = pos[:]; del pos_minus[i]
 q, p = compute_model(seq_minus, pos_minus, alphabet, w)
...
def compute_model(sequences, pos, alphabet, w):
    """
    This method compute the probability model of background and word based on data in 
    the sequences.
    """
    q = {x:[1]*w for x in alphabet}
    p = {x: 1 for x in alphabet}
    
    # Count the number of character occurrence in the particular position of word
    for i in xrange(len(sequences)):
        start_pos = pos[i]        
        for j in xrange(w):
            c = sequences[i][start_pos+j]
            q[c][j] += 1
    # Normalize the count
    for c in alphabet:
        for j in xrange(w):
            q[c][j] = q[c][j] / float( len(sequences)+len(alphabet) )
    
    # Count the number of character occurrence in background position
    # which mean everywhere except in the word position
    for i in xrange(len(sequences)):
        for j in xrange(len(sequences[i])):
            if j < pos[i] or j > pos[i]+w:
                c = sequences[i][j]
                p[c] += 1
    # Normalize the count
    total = sum(p.values())
    for c in alphabet:
        p[c] = p[c] / float(total)
    
    return q, p

Setelah mendapatkan model untuk motif dan background, kita coba hitung probabilitas untuk semua posisi motif yang mungkin pada deret i dengan cara mengalikannya dengan dua model yang sudah dihitung sebelumnya. Semua karakter dimulai dari posisi awal hingga panjang motif dikalikan dengan nilai probabilitas kemunculan karakter tersebut pada posisi yang bersangkutan.

N = len(sequences[i])
qx = [1]*(N-w+1)
px = [1]*(N-w+1)
for j in xrange(N-w+1):
    for k in xrange(w):
        c = sequences[i][j+k]
        qx[j] = qx[j] * q[c][k]
        px[j] = px[j] * p[c]

Setelah itu, kita bisa menghitung rasio antara motif dan background untuk setiap kemungkinan posisi motif pada deret tersebut. Rasio yang paling besar menandakan bahwa posisi tersebut adalah posisi yang paling mungkin sebagai awal dari motif.

Aj = [x/y for (x,y) in zip(qx, px)]
norm_c = sum(Aj)
Aj = map(lambda x: x/norm_c, Aj)

Kita lakukan sampling posisi baru berdasarkan distribusi rasio. Rasio yang lebih besar tentunya memiliki peluang untuk terpilih lebih besar. Hal ini dilakukan untuk memberikan randomness dan menghindari local maxima (walaupun mungkin kurang berhasil).

pos[i] = sample(range(N-w+1), Aj)
...
def sample(alphabet, dist):
    """ This method produce a new discrete sample list by alphabet with probability
    distribution given by dist.
    The length of alphabet and dist must be same."""
    sampl = None
    cum_dist = np.cumsum(dist)
    r = rand()
    for i in xrange(len(dist)):
        if r < cum_dist[i]:
            sampl = alphabet[i]
            break
    
    return sampl

Hal ini diulang untuk semua deret, dan dilakukan hingga tidak ada lagi perubahan posisi untuk motif.

Kode lengkapnya bisa dilihat di https://github.com/mitbal/gibbs-sampler-motif-finding pada file gibbs.py

Sekarang untuk tes algoritma kita buat satu file python baru dan masukkan kode berikut

seqs = ['muhammadiqbal', 'iqbaltawakal', 'miqbalt']
k = 5

new_pos = gibbs.sampling(seqs, k)

words = [seqs[i][new_pos[i]:new_pos[i]+k] for i in xrange(len(seqs))]
print words

Jika sukses, seharusnya program di atas mengeluarkan output seperti berikut

['iqbal', 'iqbal', 'iqbal']

Namun terkadang Gibbs sampling bisa terjebak di local maxima karena inisialisasi posisi awal yang kurang baik atau pemilihan sampling posisi baru yang kurang beruntung.

['uhamm', 'iqbal', 'iqbal']

Untuk mengurangi risiko tersebut, kita bisa menjalankan gibbs sampling beberapa kali dan menghitung berapa kali output yang dihasilkan keluar. Dengan memilih output yang paling sering keluar kita bisa meningkatkan akurasi dari Gibbs sampling. Istilahnya adalah multiple chains.

result = {}
for i in xrange(20):
    new_pos = gibbs.sampling(seqs, k)
    #print new_pos
    tnp = tuple(new_pos)    
    if tnp in result:
        result[tnp] += 1
    else:
        result[tnp] = 1

max_vote = 0
max_pos = None
for key in result:
    #print key, result[key]
    if result[key] > max_vote:
        max_pos = list(key)
        max_vote = result[key]

words = [seqs[i][max_pos[i]:max_pos[i]+k] for i in xrange(len(seqs))]
print words

Selain cara diatas, module gibbs bisa dipanggil langsung dari command line. Pertama masukkan data mengenai deret dalam sebuah file dengan format baris pertama adalah panjang motif, dan kemudian setiap barisnya adalah satu deret. Contoh file “test.txt”

3
thequickdog
browndog
dogwood

Kemudian panggil program dari terminal dengan cara

 cat test.txt | python gibbs.py

atau

 python gibbs.py test.txt

Output-nya adalah posisi motif untuk setiap deret dalam bentuk list.

Semoga berguna. Salam.

One thought on “Gibbs Sampling

Tinggalkan Balasan

Isikan data di bawah atau klik salah satu ikon untuk log in:

Logo WordPress.com

You are commenting using your WordPress.com account. Logout / Ubah )

Gambar Twitter

You are commenting using your Twitter account. Logout / Ubah )

Foto Facebook

You are commenting using your Facebook account. Logout / Ubah )

Foto Google+

You are commenting using your Google+ account. Logout / Ubah )

Connecting to %s