al2na

download.file("https://dl.dropboxusercontent.com/u/1373164/H1.chr21.chr22.rds",destfile="H1.chr21.chr22.rds",method="curl")

mbw=readRDS("H1.chr21.chr22.rds")

# it finds the optimal number of componets as 6
res=methSeg(mbw,diagnostic.plot=TRUE,maxInt=100,minSeg=10)
# however the BIC stabilizes after 4, we can also try 4 componets
res=methSeg(mbw,diagnostic.plot=TRUE,maxInt=100,minSeg=10,G=1:4)

# get segments to BED file

al2na

.Rproj.user
.Rhistory
.RData
*.Rproj
*.html

al2na

#' Read bismark coverage file as a methylKit object
#' 
#' Bismark aligner can output methylation information per base in
#' multiple different formats. This function reads coverage files,
#' which have chr,start,end, number of cytosines (methylated bases) 
#' and number of thymines (unmethylated bases).
#' 
#' @param location a list or vector of file paths to coverage files
#'     
#' @param sample.id a list or vector of sample ids

al2na

library(data.table)
library(genomation)
library(ggplot2)


#' Annotate given ranges with genomic features
#' 
#' The function annotates a target GRangesList or GRanges object as overlapping
#' or not with 
#' the elements of named GRangesList. This is useful to annotate your regions

al2na

p300.url="http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsHaibH1hescP300V0416102UniPk.narrowPeak.gz"

Nanog.url="http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsHaibH1hescNanogsc33759V0416102UniPk.narrowPeak.gz"

SP1.url="http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsHaibH1hescSp1Pcr1xUniPk.narrowPeak.gz"

chrHMM.url="http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmH1hescHMM.bed.gz"
#http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgSegmentation/wgEncodeAwgSegmentationChromhmmH1hesc.bed.gz"

# get chromHMM annotation and make a list out of it