One of the features of methpipe is using a Hidden Markov Model to determine if sites are hypo-methylatyed or not. I figured I’d give that a whirl with the Day 10 data, since it was finished.
First, according to the methpipe documentation we want to restrict our data to only the CpG sites. They offer a handy command “symmetric-cpgs” to do this.
Set working directory, and get a list of the .meth files located in it.
setwd("~/Documents/methpipe-Day10-countfiles")
meth.files <- list.files(pattern = "*.meth")
meth.files [1] "EPI-103.meth" "EPI-104.meth" "EPI-111.meth" "EPI-113.meth" "EPI-119.meth" "EPI-120.meth" "EPI-127.meth" "EPI-128.meth" "EPI-135.meth"
[10] "EPI-136.meth" "EPI-143.meth" "EPI-145.meth"Run the symmetric-cpgs command.
for(i in 1:length(meth.files)) {
system(paste0("/home/shared/methpipe/methpipe-3.4.3/bin/symmetric-cpgs -o ", substr(meth.files[i], 1, 7), "-CpG.meth ", meth.files[i]))
}cpg.files <- list.files(pattern = "*-CpG.meth")
cpg.files [1] "EPI-103-CpG.meth" "EPI-104-CpG.meth" "EPI-111-CpG.meth" "EPI-113-CpG.meth" "EPI-119-CpG.meth" "EPI-120-CpG.meth" "EPI-127-CpG.meth"
[8] "EPI-128-CpG.meth" "EPI-135-CpG.meth" "EPI-136-CpG.meth" "EPI-143-CpG.meth" "EPI-145-CpG.meth"Next, run the hmr program, which identifies hypomethylated regions using a hidden markov model by loci.
for(i in 1:length(meth.files)) {
system(paste0("/home/shared/methpipe/methpipe-3.4.3/bin/hmr -o ", substr(cpg.files[i], 1, 7), "-CpG.hmr ", cpg.files[i]))
}Make sure the HMR files exist, and look at a few of them.
hmr.files <- list.files(pattern = "*.hmr")
hmr.files [1] "EPI-103-CpG.hmr" "EPI-104-CpG.hmr" "EPI-111-CpG.hmr" "EPI-113-CpG.hmr" "EPI-119-CpG.hmr" "EPI-120-CpG.hmr" "EPI-127-CpG.hmr"
[8] "EPI-128-CpG.hmr" "EPI-135-CpG.hmr" "EPI-136-CpG.hmr" "EPI-143-CpG.hmr" "EPI-145-CpG.hmr"system(paste0("head ", hmr.files[1]))scaffold10000 244 3634 HYPO0 45 +
scaffold100042 1235 9910 HYPO1 92 +
scaffold100043 518 3011 HYPO2 73 +
scaffold100043 5922 13424 HYPO3 117 +
scaffold10005 171 5692 HYPO4 65 +
scaffold10005 6736 18646 HYPO5 141 +
scaffold10005 22391 39276 HYPO6 235 +
scaffold100063 361 8039 HYPO7 92 +
scaffold100077 526 7250 HYPO8 109 +
scaffold100105 589 4022 HYPO9 65 +system(paste0("head ", hmr.files[2]))scaffold10000 565 4240 HYPO0 42 +
scaffold10000 5463 8904 HYPO1 39 +
scaffold100042 1343 9910 HYPO2 88 +
scaffold100043 518 3578 HYPO3 84 +
scaffold100043 5943 12660 HYPO4 103 +
scaffold10005 313 21174 HYPO5 253 +
scaffold10005 22441 39091 HYPO6 269 +
scaffold100077 363 7250 HYPO7 123 +
scaffold100105 589 4154 HYPO8 65 +
scaffold100105 5433 10596 HYPO9 36 +Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).