Dinh/Dinh 2012/NOTES/2012-1-26

From ZhangLabWiki
Jump to navigation Jump to search

HMM to determine UMR, LMR, and HMR[edit]

Theory of HMM[edit]

  • A Markov model has:
    • observable states: 1,2,...,N
    • observed sequence: q1,q2,...,qr
    • first order Markov assumption: the next state depends only on the current state and not on the sequence of events that preceded it.
    • at each state, there are emission probabilities of emission any one of the observed symbols
    • at each state, there are transition probabilities of remaining in the same state or changing to another state.
  • A Hidden Markov model has:
    • states are not observable, and so the observed sequence gives the probability of states.
    • give N number of hidden states, a sequence of length M, and transition and emission probability matrices, find the most likely sequence of states.
      • this problem is solved using a viterbi path algorithm (dynamic programming).

Applied HMM[edit]

  • install R package "RHmm"
  • Let the hiddden states be UMR, LMR, and HMR for a sequence of CpG sites
  • Let the observed sequence be the fractional methylation level (0-1)
  • Let the emission probability be a normal distribution with mean and standard deviation for each state.
  • First, to build the HMM, the Baum-Welch algorithm is used to estimate the transition and emission probabilities
myFit = HMMFit(obs=myObs, nStates=3, control=list(verbose=1, init="KMEANS"), asymptCov=TRUE, asymptMethod='optim')
    • To better build the model, I will use estimates from the entire chromosome 20 and from multiple tissue types.
sort -k2,2n /home/nplongth/Noi_scratch/HMM-LMRs_2012_01_25/wgbs.BED/h1.wgbs.chro20.BED.txt | awk '{print $2"\t"$4}' > h1.seq
(repeat for other tissue types, then remove the first header line for every *.seq file).
R --max-ppsize=500000
    • Then run:
>source("calculateHMM.R")
./combineStatesandPos.pl ads.seq ads.states | sed 's/S1/HMR/g' | sed 's/S2/LMR/g' | sed 's/S3/UMR/g' | sed 's/C1/210,0,0/g' | sed 's/C2/120,120,0/g' | sed 's/C3/0,240,0/g' > ads.combo.BED
Result from genome browser:
track: ads.seq = HMM track, with bright green (UMR), dark green (LMR), and red (HMR)
track: Sample = CpG methylation level green (0) to red (1)
File:Hmm ads.png