Whatsapp Chat Analyzer

Jadi ceritanya mahasiswa asal Indonesia yang tinggal di Stockholm menggunakan whatsapp sebagai media komunikasi untuk saling berbagi informasi hingga keluh kesah. Karena jumlah pesertanya yang memang tidak terlalu banyak (20+++), saya rasa, dibanding grup Facebook atau mailing list, group chat whatsapp memang lebih cocok untuk komunikasi yang membutuhkan respons yang relatif cepat.

Berawal dari rasa ingin tahu akan esensi dari pembicaraan di grup whatsapp tersebut, saya pun penasaran dan ingin menganalisis kontennya. Bak gayung bersambut, ternyata whatsapp menyediakan fitur untuk mengarsip pembicaraan dan mengirimkannya lewat email. Caranya cukup tahan percakapan yang ingin diarsip. Kemudian pilih ‘Email conversation’ > ‘Without media’. Gambar di bawah menunjukkan langkah-langkah tersebut. Nama dan tulisan disamarkan untuk menjaga privasi.

archive conversationwithout media

Hasilnya adalah sebuah file text. Setiap pesan chat ditulis pada satu baris dengan format “Tanggal – User: Pesan” (tanpa tanda kutip). Kode python dibawah digunakan untuk membaca file tersebut.

def read(filename):
    """ Read text file containing whatsapp chat and return the list of list of time, author, and its text
    :param filename: the filename of the chat text file
    :return: chat 2d list
    """
    chat = []
    with open(filename, 'r') as f:
        for line in f:
            lines = line.split(' - ')  # Divide between date and the rest
            if len(lines) > 1:
                lines2 = lines[1].split(': ')  # Divide between user and text
                if len(lines2) > 1:
                    speaker = lines2[0]
                    text = lines2[1]
                else:
                    speaker = ''
                    text = lines2[0]
                timestamp = lines[0]
            else:
                timestamp = ''
                speaker = ''
                text = lines[0]
            chat += [[timestamp, speaker, text]]
    return chat

Sip, setelah kita punya datanya, kita bisa geledah dan obrak-abrik isinya untuk dianalisis lebih lanjut.

Statistika dasar

Pertama-pertama. kita ingin tahu ada berapa banyak baris dalam percakapan.

date

Hmmm, ternyata ada sekitar 137 pesan per harinya. Lumayan aktif juga.

Kemudian kita ingin tahu siapa sih yang paling sering bunyi dan siapa yang paling senyap. Kita hitung total berapa banyak pesan yang ditulis per orangnya.

user frequency
Nama dan nomor telepon disamarkan untuk menjaga privasi

Wah, ternyata ada perbedaan yang cukup jomplang antara yang paling sering mengirim pesan dengan yang paling tidak sering mengirim pesan. Bahkan yang nomor 1 jumlah pesannya sekitar 2 kali dari yang nomor 2. Distribusi ini sepertinya mengikuti Zipf’s Law.

Frekuensi kata

Fitur utama dalam analisis dokumen biasanya adalah frekuensi kata. Model ini biasanya dikenal dengan nama Bag-of-Words. Kalau misalkan fitur ini dirasa kurang cukup ekspresif, biasanya ditambah lagi menjadi n-grams model, dimana kata dihitung kemunculannya bersama kata lain.

Dalam analisis dokumen juga biasanya terdapat kata-kata yang frekuensinya tinggi namun tidak memberi makna lebih kepada tulisan, seperti kata hubung dan teman-temannya. Kumpulan kata-kata tersebut sering disebut dengan istilah stopwords. Stopwords ini biasanya dihilangkan sebelum analisis agar proses lebih ringkas dan akurat. Saya memakai daftar kata stopwords dari https://sites.google.com/site/kevinbouge/stopwords-lists yang ternyata mengambil dari thesis “A Study of Stemming Effects on Information Retrieval in Bahasa Indonesia“.

Sekarang kita coba hitung frekuensi kata dalam chat minus stopwords. Histogram hasilnya bisa dilihat di bawah.

word frequency

Hmm, menarik. Ternyata kata-kata yang paling banyak keluar adalah kata sambung yang salah ketik seperti ‘yg‘ dan ‘d‘.  Kata-kata lain yang banyak muncul adalah nama panggilan ke penghuni grup yang lain. Selain itu, ternyata grup ini cukup suka tertawa, terlihat dari 3 kemunculan ekspresi tawa yaitu ‘hahaha‘, ‘wkwk‘, dan ‘haha‘. Selain itu, ternyata grup ini rajin olah raga juga, dilihat dari kata ‘badminton‘ yang muncul lebih dari 200 kali.

Emoji

if a picture paints a thousand words then why can’t I paint you“. Kurang lebih itulah sepenggal lagu If dari Bread. Penggunaan emoji atau emoticon sangat kentara dalam chatting karena satu gambarnya bisa mengekspresikan hal yang sulit jika dicoba disampaikan dengan kata-kata. Penggunaan emoji ini sangat pesat hingga menimbulkan sentimen kalau nanti kita bisa kembali ke zaman Mesir kuno dengan tulisan hieroglyph-nya. Karena umumnya penggunaan emoji ini, saya pun mencoba menghitung frekuensi dari tiap-tiap emoji.

Yang agak menyulitkan mungkin adalah emoji ini direpresentasikan dalam unicode dan di-encode dalam UTF-8 sehingga harus dibuat fungsi konversinya dulu antar keduanya. Gambar untuk plot didapat dari https://github.com/github/gemoji/ yang saya dapat link-nya dari apps.timwhitlock.info/emoji/tables/unicode. Kesulitan lain adalah ternyata Python Image Library punya sedikit untuk png dengan transparency. Apa boleh buat, gambar emoji-nya pun harus saya konversi ke jpg dengan utilitas mogrify dari imagemagick. Histogram hasilnya bisa dilihat di gambar di bawah.

emoji frequency

Hmm, menarik. Bisa dilihat, emoji teratas adalah senyum sambil berkeringat yang saya baca sebagai ekspresi dari ‘hehe’. Sisanya adalah berbagai ekspresi untuk senyum dan tertawa. Tetapi emoji ke 3 yang agak aneh. Mata. Ada apa dengan mata? Emoji ke 8, ‘see no evil monkey’ juga aneh. Apakah ada hubungannya dengan emoji ke 3?

Topik

Target berikutnya adalah Topic modelling. Lazimnya, setahu saya ini memakai teknik Latent Dirichlet Allocation (LDA). Karena saya belum pernah implementasi teknik ini sebelumnya dan kekurangan waktu untuk mencoba, saya mencoba mencari apakah ada library lain siap pakai. Saya ketemu situs startup untuk klasifikasi teks http://www.monkeylearn.com/. Di situs tersebut terdapat servis untuk melakukan klasifikasi topik dokumen secara general. Untuk free account diberikan 1000 kali panggilan API. Yah, cukuplah untuk proyek kecil-kecilan ini.

monkeylearn1monkeylearn2

Karena model-nya dibangun dari korpus Bahasa Inggris, maka teks harus diterjemahkan dulu ke bahasa tersebut. Hal ini sebenarnya bisa mengurangi akurasi, akan tetapi akurasi bukan tujuan utama yang dicari. Pilihan pertama saya, Google Translate API, ternyata tidak gratis. Ini tentunya cukup memberatkan saya yang cuma mahasiswa. Saya pun memakai layanan alternatif dari Yandex, search engine buatan Rusia.

topic

Dan inilah hasilnya. Ternyata percakapan di grup whatsapp ini dikategorikan ke ‘Entertainment & Recreation’. Hmm, bisa jadi. Yang agak aneh adalah kategori kedua, ‘Anime’, yang merupakan subkategori dari Entertainment di atas karena seingat saya tidak ada diskusi tentang anime sama sekali di grup. Hmm, inilah sulitnya memakai model buatan orang. Sulit untuk dianalisis.

Sebenarnya masih ada lagi yang mau saya coba, seperti menggunakan t-SNE untuk melihat kemiripan antar user. Namun, apa daya TTM (thesis telah memanggil).

Kode lengkap dapat dilihat di https://github.com/mitbal/wca. Cukup ikuti petunjuk di README untuk menjalankan program.

Semoga bermanfaat. Salam.

Baseline Wander Removal dengan Wavelet

Pada post kali ini, kita akan melakukan proses menghilangkan baseline wander (baseline wander removal) pada sinyal EKG dengan menggunakan transformasi wavelet seperti yang ditulis oleh Sargolzaei et al pada paper dengan judul “A new robust wavelet based algorithm for baseline wandering cancellation in ECG signals“.

Latar belakang
Elektrokardiogram (EKG) adalah pembacaan sinyal elektrik jantung dengan cara menempelkan elektroda ke posisi tertentu pada tubuh dan kemudian membaca nilai perbedaan potensial listrik alias tegangan yang dihasilkan. Dari pembacaan ini kita bisa mencari tahu mengenai kondisi jantung. Akan tetapi, sinyal yang baru dibaca ini tidak lepas dari noise sehingga bisa mengganggu proses diagnosis. Sumber noise ini berasal dari berbagai macam, seperti elektroda dari elektroda itu sendiri atau jika tubuh bergerak ketika dilakukan pengukuran. Salah satu kondisi noise yang sering menyerang adalah situasi yang disebut baseline wander. Baseline wander terjadi apabila sinyal EKG tidak lurus pada sumbu x, malah naik turun. Contoh bisa dilihat di gambar 1.

bwr
Contoh sinyal yang terkena baseline wander dan hasilnya setelah dihilangkan. Gambar diambil tanpa izin dari paper yang dirujuk di atas.

Baseline wander removal (BWR) adalah salah satu tahap preprocessing pada sinyal EKG untuk menghilangkan baseline drift ini. Ada beberapa teknik yang bisa dipakai, seperti melakukan filtering. Pada post ini, kita akan mengimplementasikan salah satu teknik yang memanfaatkan transformasi wavelet. Penjelasan dasar mengenai wavelet salah satunya bisa dilihat di http://users.rowan.edu/~polikar/WAVELETS/WTtutorial.html. Kelebihan dari teknik ini adalah kita tidak perlu menentukan parameter seperti ketika menggunakan high pass filter (frequency cut-off) sehingga metode kita bekerja secara non-supervised.

Sinyal EKG yang sudah bersih dari noise semacam ini kemudian bisa digunakan untuk berbagai macam kegunaan. Untuk pembahasan mengenai kegunaan sinyal EKG untuk mengenali tipe detak arrhythmia atau diagnosis tahapan tidur, bisa merujuk ke paper ini dan ini.

Algoritma
Berikut adalah tahapan dalam algoritma ini. Pertama, kita lakukan dekomposisi kepada sinyal asli dengan transformasi wavelet. Dalam hal ini kita memilih Daubechies orde 4 sebagai fungsi basisnya. Sinyal didekomposisi menjadi bagian frekuensi rendah/aproksimasi dan frekuensi tinggi/detail. Kemudian kita hitung nilai energi pada sinyal frekuensi tinggi. Kita cari kondisi dimana nilai energi pada level dekomposisi tersebut lebih rendah dari pada nilai pada level dekomposisi sebelumnya dan sesudahnya (alias lokal minima). Setelah kita temukan level tersebut, kita rekonstruksi sinyal aproksimasi dari level ini dengan membuang nilai pada sinyal frekuensi tinggi (atau dijadikan 0 semua). Sinyal hasil rekonstruksi ini kita sebut dengan baseline. Untuk menghilangkan baseline wander pada sinyal asli, maka kita kurangkan sinyal asli dengan sinyal baseline.

level
Gambar atas: plot nilai energi pada sinyal detail pada berbagai level dekomposisi. Tanda panah menunjukkan level ketika nilainya adalah lokal minima. Gambar bawah: Sinyal asli, baseline, dan sinyal asli yang sudah dikurang baseline. Lagi-lagi diambil dari paper rujukan di atas.

 

Implementasi
Pertama, kita implementasi metode konvolusi antara dua sinyal.

def conv(x, h):
    """ Perform the convolution operation between two input signals. The output signal length
    is the sum of the lenght of both input signal minus 1."""
    length = len(x) + len(h) - 1
    y = [0]*length

    for i in xrange(len(y)):
        for j in xrange(len(h)):
            if i-j >= 0 and i-j < len(x):
                y[i] += h[j] * x[i-j]

    return y

Lalu kita implementasi metode dekomposisi wavelet dan jangan lupa deklarasi koefisien basis fungsi yang dipakai. Dalam kasus ini adalah Daubechies dengan 4 koefisien. Dekomposisi dilakukan dengan cara mengkonvolusi sinyal asli dengan koefisien low pass and high pass. Sinyal keluaran dari proses ini kemudian di-downsampling menjadi setengahnya dengan cara hanya mengambil nilai pada posisi ganjil atau genap saja. Hal ini dilakukan berulang terhadap sinyal keluaran low pass (atau disebut aproksimasi) sebanyak parameter level.

c0 = (1+sqrt(3))/(4*sqrt(2))
c1 = (3+sqrt(3))/(4*sqrt(2))
c2 = (3-sqrt(3))/(4*sqrt(2))
c3 = (1-sqrt(3))/(4*sqrt(2))

def db4_dec(x, level):
""" Perform the wavelet decomposition to signal x with Daubechies order 4 basis function as many as specified level"""

    # Decomposition coefficient for low pass and high pass
    lpk = [c0, c1, c2, c3]
    hpk = [c3, -c2, c1, -c0]

    result = [[]]*(level+1)
    x_temp = x[:]
    for i in xrange(level):
        lp = conv(x_temp, lpk)
        hp = conv(x_temp, hpk)

        # Downsample both output by half
        lp_ds=[0]*(len(lp)/2)
        hp_ds=[0]*(len(hp)/2)
        for j in xrange(len(lp_ds)):
            lp_ds[j] = lp[2*j+1]
            hp_ds[j] = hp[2*j+1]

        result[level-i] = hp_ds
        x_temp = lp_ds[:]

    result[0] = lp_ds
    return result

Fungsi rekontruksi digunakan untuk mencari baseline dari sinyal asal. Rekonstruksi bekerja dengan melakukan konvolusi dengan koefisien rekonstruksi pada kedua sinyal low pass dan high pass, melakukan upsampling dengan menyelipkan 0 di setiap nilai pada sinyal dan kemudian masing-masing posisi pada sinyal saling dijumlahkan.

def db4_rec(signals, level):
    """ Perform reconstruction from a set of decomposed low pass and high pass signals as deep as specified level"""

    # Reconstruction coefficient
    lpk = [c3, c2, c1, c0]
    hpk = [-c0, c1, -c2, c3]

    cp_sig = signals[:]
    for i in xrange(level):
        lp = cp_sig[0]
        hp = cp_sig[1]

        # Verify new length
        length = 0
        if len(lp) > len(hp):
            length = 2*len(hp)
        else:
            length = 2*len(lp)

        # Upsampling by 2
        lpu = [0]*(length+1)
        hpu = [0]*(length+1)
        index = 0
        for j in xrange(length+1):
            if j%2 != 0:
                lpu[j] = lp[index]
                hpu[j] = hp[index]
                index += 1

        # Convolve with reconstruction coefficient
        lpc = conv(lpu, lpk)
        hpc = conv(hpu, hpk)

        # Truncate the convolved output by the length of filter kernel minus 1 at both end of the signal
        lpt = lpc[3:-3]
        hpt = hpc[3:-3]

        # Add both signals
        org = [0]*len(lpt)
        for j in xrange(len(org)):
            org[j] = lpt[j] + hpt[j]

        if len(cp_sig) > 2:
            cp_sig = [org]+cp_sig[2:]
        else:
            cp_sig = [org]

    return cp_sig[0]

Method calcEnergy menghitung nilai energi dari sebuah sinyal berdasarkan definisinya yaitu jumlahan kuadrat sinyal di tiap titik

def calcEnergy(x):
    """ Calculate the energy of a signal which is the sum of square of each points in the signal."""
    total = 0
    for i in x:
        total += i*i
    return total

Kemudian method bwr berikut adalah implementasi dari algoritma yang dijabarkan di atas.

def bwr(raw):
    """ Perform the baseline wander removal process against signal raw. The output of this method is signal with correct baseline
    and its baseline """
    en0 = 0
    en1 = 0
    en2 = 0
    n = 0

    curlp = raw[:]
    num_dec = 0
    last_lp = []
    while True:
        print 'Iterasi ke' + str(num_dec+1)
        print len(curlp)

        # Decompose 1 level
        [lp, hp] = db4_dec(curlp,1)

        # Shift and calculate the energy of detail/high pass coefficient
        en0 = en1
        en1 = en2
        en2 = calcEnergy(hp)
        print en2

        # Check if we are in the local minimum of energy function of high-pass signal
        if en0 > en1 and en1 < en2:
            last_lp = curlp
            break

        curlp = lp[:]
        num_dec = num_dec+1

    # Reconstruct the baseline from this level low pass signal up to the original length
    base = last_lp[:]
    for i in xrange(num_dec):
        base = db4_rec([base,[0]*len(base)], 1)

    # Correct the original signal by subtract it with its baseline
    ecg_out = [0]*len(raw)
    for i in xrange(len(raw)):
        ecg_out[i] =  raw[i] - base[i]

    return (base, ecg_out)

Contoh
Untuk contoh kita ambil data dari situs Physionet khususnya MIT-BIH Arrhythmia database. Data diambil dari salah satu record pasien dengan kode nomor 101. Sinyal EKG diambil sepanjang 1 menit. Sinyal EKG yang diambil berasal dari lead II dan V5.

Contoh kode untuk memanggil modul dan fungsi yang sudah dibuat di atas adalah sebagai berikut

import bwr
import matplotlib.pyplot as plt

# Read input csv file from physionet
f = open('samples1.csv', 'r')
lines = f.readlines()
f.close()

# Discard the first two lines because of header. Takes either column 1 or 2 from each lines (different signal lead)
signal = [0]*(len(lines)-2)
for i in xrange(len(signal)):
	signal[i] = float(lines[i+2].split(',')[1])

# Call the BWR method
(baseline, ecg_out) = bwr.bwr(signal)

plt.subplot(2,1,1)
plt.plot(signal, 'b-')
plt.plot(baseline, 'r-')

plt.subplot(2,1,2)
plt.plot(ecg_out, 'b-')
plt.show()

Berikut contoh pertama. Garis merah pada plot merupakan baseline dari sinyal EKG.

samples1

Berikut adalah contoh kedua.

samples2

Contoh berikut diambil dari lead V5.

samples3

Kode lengkap dan sampel sinyal dapat dilihat di: https://github.com/mitbal/py-bwr

Semoga berguna. Salam.

UPGMA

Pada post kali ini kita akan mengimplementasi algoritma UPGMA atau Unweighted Pair Group Method with Arithmetic Mean untuk melakukan clustering.

Algoritma
UPGMA bekerja dengan prinsip yang sederhana. Pada setiap iterasi, pilih pasangan point dengan point, atau point dengan cluster, atau cluster dengan cluster, yang memiliki jarak yang terpendek. Gabung kedua pasangan ini kedalam satu cluster. Hal ini dilakukan terus menerus hingga jumlah cluster berkurang menjadi jumlah yang diinginkan. Untuk menghitung jarak antar datapoint, kita bisa menggunakan berbagai kriteria jarak. Pada implementasi ini, kita menggunakan euclidean distance untuk menghitung jarak. Sedangkan untuk menghitung jarak antar cluster, kita hitung dengan cara merata-ratakan jarak antar semua pasang point dari cluster pertama ke cluster kedua, atau dengan kata lain:

\frac{1}{|A|.|B|} \sum_{x \in A} \sum_{y \in B} d(x,y)

Teknik UPGMA ini bisa digolongkan kepada teknik hierarchical clustering karena kita bisa melihat hirarki sebuah cluster yang dibentuk dari cluster-cluster lain yang lebih kecil.

Implementasi
Kita implementasikan UPGMA dalam bahasa Python. Pertama, kita butuh struktur data untuk menyimpan point yang terdapat pada sebuah cluster.

class Node:
    def __init__(self, p):
        self.points = p
        self.right = None
        self.left = None

Fungsi UPGMA menerima dua parameter. Parameter pertama adalah set yang berisi datapoints sebanyak n dengan masing-masing point memiliki dimensi d. Parameter kedua adalah jumlah cluster yang kita inginkan.

def upgma(points, k):
    """ Cluster based on distance matrix dist using Unweighted Pair Group Method with Arithmetic Mean algorithm up to k cluster"""

    # Initialize each cluster with one point
    nodes = []
    n = len(points)
    for i in xrange(n):
        node = Node([points[i]])
        nodes = nodes + [node]

    # Iterate until the number of clusters is k
    nc = n
    while nc > k:
        # Calculate the pairwise distance of each cluster, while searching for pair with least distance
        c1 = 0; c2 = 0; i1 = 0; i2 = 0;
        sdis = 9999999999
        for i in xrange(nc):
            for j in xrange(i+1, nc):
                dis = euclidistance(nodes[i], nodes[j])
                if dis < sdis:
                    sdis = dis
                    c1 = nodes[i]; c2 = nodes[j];
                    i1 = i; i2 = j;
        # Merge these two nodes into one new node
        node = Node(c1.points + c2.points)
        node.left = c1; node.right = c2;
        
        #Remove the previous nodes, and add the new node
        new_nodes = []
        for i in xrange(nc):
            if i != i1 and i != i2:
                new_nodes = new_nodes + [nodes[i]]
        new_nodes = new_nodes + [node]
        nodes = new_nodes[:]
        nc = nc - 1

    return nodes

Kemudian terakhir, kita definisikan fungsi jarak yang kita pakai. Selain jarak euclidean, kita bisa yang digunakan kriteria jarak lainnya seperti Manhattan distance atau Chebyshev distance.

import math

def euclidistance(c1, c2):
    """ Calculate the distance between two cluster """
    dist = .0
    n1 = len(c1.points)
    n2 = len(c2.points)
    for i in xrange(n1):
        for j in xrange(n2):
            p1 = c1.points[i]
            p2 = c2.points[j]
            dim = len(p1)
            d = 0
            for k in xrange(dim):
                d = d + (p1[k]-p2[k])**2
            d = math.sqrt(d)
            dist = dist + d
    dist = dist / (n1*n2)
    return dist

Contoh
Untuk contoh pertama kita pakai dataset sintesis yang dihasilkan oleh dua buah distribusi normal multidimensi.

import upgma
import random
import matplotlib.pyplot as plt
import math

# Example 1
datapoints = [(random.normalvariate(2.5, 1.0), random.normalvariate(1.5,1.0)) for i in xrange(100)] + \
				[(random.normalvariate(-1, 0.5), random.normalvariate(3,0.5)) for i in xrange(100)]

# Plot datapoints before clustering
plt.plot([x for (x,y) in datapoints], [y for (x,y) in datapoints], 'k^')
plt.show()

# Cluster the data
nodes = upgma.upgma(datapoints, 2)
plt.plot([x[0] for x in nodes[0].points], [x[1] for x in nodes[0].points], 'b*')
plt.plot([x[0] for x in nodes[1].points], [x[1] for x in nodes[1].points], 'ro')
plt.show()

dan berikut keluaran hasil clustering jika kita pilih 2 sebagai jumlah cluster.
nocluster1cluster1

Contoh kedua kita ambil dataset old faithful geyser (http://www.stat.cmu.edu/~larry/all-of-statistics/=data/faithful.dat).

f = open('faithful.dat.txt', 'r')
lines = f.readlines()
f.close()

datapoints = []
for line in lines:
	tokens = line.split()
	datapoints += [(float(tokens[1]), float(tokens[2]))]

Karena fitur kedua berbeda skalanya dengan fitur pertama (puluhan berbanding dengan satuan), kita lakukan normalisasi z-score pada tiap masing-masing fitur dengan cara mengurangi dengan mean distribusi dan kemudian dibagi dengan standar deviasi.

avg1 = sum([x for (x,y) in datapoints])
avg2 = sum([y for (x,y) in datapoints])
centered_datapoints = map(lambda (x,y): (x-avg1, y-avg2), datapoints)
std1 = math.sqrt(sum(map(lambda x: x*x, [x for (x,y) in centered_datapoints])))
std2 = math.sqrt(sum(map(lambda x: x*x, [y for (x,y) in centered_datapoints])))
normalized_datapoints = map(lambda (x,y): (x/std1, y/std2), centered_datapoints)

Hasil clustering-nya adalah sebagai berikut.

# Before clustering
plt.plot([x for (x,y) in normalized_datapoints], [y for (x,y) in normalized_datapoints], 'k^')
plt.show()

# Cluster the data
nodes = upgma.upgma(normalized_datapoints, 2)
plt.plot([x[0] for x in nodes[0].points], [x[1] for x in nodes[0].points], 'b*')
plt.plot([x[0] for x in nodes[1].points], [x[1] for x in nodes[1].points], 'ro')
plt.show()

nocluster2cluster2

Kode lengkap dapat dilihat di: https://github.com/mitbal/py-upgma

Semoga berguna. Salam.

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.