Dinh/Dinh 2012/NOTES/2012-1-26
Jump to navigation
Jump to search
HMM to determine UMR, LMR, and HMR[edit]
- Idea from http://www.nature.com/nature/journal/v480/n7378/extref/nature10716-s1.pdf
- (note link to supplementary of the paper only!)
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 script to run HMM: File:CalculateHMM.R.txt
- First start R with increased ppsize to avoid stack overflow error.
R --max-ppsize=500000
- Then run:
>source("calculateHMM.R")
- Outputs from R: File:CalculateHMM.Rout.txt]]
- Script to create BED files: File:CombineStatesandPos.txt
./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