### R code from vignette source 'vignettes/biomvRCNS/inst/doc/biomvRCNS.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: loadlib
###################################################
library(biomvRCNS)


###################################################
### code chunk number 2: setwidth
###################################################
options(width = 95)


###################################################
### code chunk number 3: coriellGR
###################################################
data('coriell', package='biomvRCNS')
head(coriell, n=3)
xgr<-GRanges(seqnames=paste('chr', coriell[,2], sep=''), 
	IRanges(start=coriell[,3], width=1, names=coriell[,1]))
values(xgr)<-DataFrame(coriell[,4:5], row.names=NULL)
xgr<-sort(xgr)
head(xgr, n=3)


###################################################
### code chunk number 4: coriellHsmm
###################################################
rhsmm<-biomvRhsmm(x=xgr, maxbp=4E4, J=3, soj.type='gamma', 
	emis.type='norm', prior.m='quantile', grp=c(1,2))


###################################################
### code chunk number 5: coriellHsmmres
###################################################
show(rhsmm)


###################################################
### code chunk number 6: coriellHsmmplot
###################################################
biomvRGviz(exprgr=xgr[seqnames(xgr)=='chr1', 'Coriell.13330'], 
	seggr=rhsmm@res[mcols(rhsmm@res)[,'SAMPLE']=='Coriell.13330'], 
	regionID='region',  tofile=FALSE)


###################################################
### code chunk number 7: coriellSeg
###################################################
rseg<-biomvRseg(x=xgr, maxbp=4E4, maxseg=10, family='norm', grp=c(1,2))


###################################################
### code chunk number 8: coriellSegres
###################################################
head(rseg@res)


###################################################
### code chunk number 9: coriellMGMR
###################################################
rmgmrh<-biomvRmgmr(xgr, q=0.9, high=T, maxgap=1000, minrun=2500, grp=c(1,2))
rmgmrl<-biomvRmgmr(xgr, q=0.1, high=F, maxgap=1000, minrun=2500, grp=c(1,2))
res<-c(rmgmrh@res, rmgmrl@res)


###################################################
### code chunk number 10: buildENCODEcgr (eval = FALSE)
###################################################
## winsize<-25
## cgr<-GRanges("chr17", strand='-', 
## 	IRanges(start=seq(7560001, 7610000, winsize), width =winsize))
## bf<-system.file("extdata", "encodeFiles.txt", package = "biomvRCNS")
## bamfiles<-read.table(bf, header=T, stringsAsFactors=F)
## library(Rsamtools)
## which<-GRanges("chr17", IRanges(7560001, 7610000))
## param<-ScanBamParam(which=which, what=scanBamWhat())
## for(i in seq_len(nrow(bamfiles))){
## 	frd<-scanBam(bamfiles[i,1], param=param)
## 	frdgr<-GRanges("chr17", strand=frd[[1]]$strand,
## 		IRanges(start=frd[[1]]$pos , end = frd[[1]]$pos+frd[[1]]$qwidth-1))
## 	mcols(cgr)<-DataFrame(mcols(cgr), DOC=countOverlaps(cgr, frdgr))
## }


###################################################
### code chunk number 11: buildENCODEgmgr (eval = FALSE)
###################################################
## af<-system.file("extdata", "gmodTP53.csv", package = "biomvRCNS")
## gtfsub<-read.table(af, fill=T, stringsAsFactors=F)
## idx<-gtfsub[,3]=='CDS' | gtfsub[,3]=='UTR'
## gmgr<-GRanges("chr17", IRanges(start=gtfsub[idx, 4], end=gtfsub[idx, 5], 
## 	names=gtfsub[idx, 13]), strand='-', TYPE=gtfsub[idx, 3])


###################################################
### code chunk number 12: poolENCODEcgr
###################################################
data(encodeTP53)
cgr<-encodeTP53$cgr
gmgr<-encodeTP53$gmgr
mcols(cgr)<-DataFrame(
	Gm12878=1+rowSums(as.matrix(mcols(cgr)[,1:2])), 
	K562=1+rowSums(as.matrix(mcols(cgr)[,3:4])) )


###################################################
### code chunk number 13: ENCODEHsmmTxDbSojourn
###################################################
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb<-TxDb.Hsapiens.UCSC.hg19.knownGene	
rhsmm<-biomvRhsmm(x=cgr, xAnno=txdb, maxbp=1E3, soj.type='gamma', 
	emis.type='pois', prior.m='quantile', q.alpha=0.01)


###################################################
### code chunk number 14: showENCODEHsmm
###################################################
rhsmm@res[mcols(rhsmm@res)[,'STATE']=='exon']


###################################################
### code chunk number 15: plotENCODEHsmmG
###################################################
g<-mcols(rhsmm@res)[,'STATE']=='exon' & mcols(rhsmm@res)[,'SAMPLE']=='Gm12878'
biomvRGviz(exprgr=cgr[,'Gm12878'], gmgr=gmgr, 
	seggr=rhsmm@res[g], plotstrand='-', regionID='TP53',  tofile=FALSE)


###################################################
### code chunk number 16: plotENCODEHsmmK
###################################################
k<-mcols(rhsmm@res)[,'STATE']=='exon' & mcols(rhsmm@res)[,'SAMPLE']=='K562'
biomvRGviz(exprgr=cgr[,'K562'], gmgr=gmgr, 
  seggr=rhsmm@res[k], plotstrand='-', regionID='TP53', tofile=FALSE)


###################################################
### code chunk number 17: findnew
###################################################
nK2gm<-findOverlaps(rhsmm@res[k], gmgr)@queryHits
nK2G<-findOverlaps(rhsmm@res[k], rhsmm@res[g])@queryHits
rhsmm@res[k][setdiff(seq_len(sum(k)), unique(c(nK2G, nK2gm)))]


###################################################
### code chunk number 18: ENCODEothers
###################################################
rseg<-biomvRseg(x=cgr, maxbp=1E3, maxseg=20, family='pois')
rmgmr<-biomvRmgmr(x=cgr, q=0.99, maxgap=50, minrun=100)


###################################################
### code chunk number 19: variodata
###################################################
data(variosm)
head(variosm, n=3)


###################################################
### code chunk number 20: varioHsmmrun
###################################################
rhsmm<-biomvRhsmm(x=variosm, maxbp=100, prior.m='quantile', r.var=0.05)


###################################################
### code chunk number 21: finddmr
###################################################
hiDiffgr<-rhsmm@res[mcols(rhsmm@res)[,'STATE']!=2 
	& mcols(rhsmm@res)[,'SAMPLE']=='meth.diff']
loPgr<-rhsmm@res[mcols(rhsmm@res)[,'STATE']==1
	& mcols(rhsmm@res)[,'SAMPLE']=='p.val']
DMRs<-intersect(hiDiffgr, loPgr)
idx<-findOverlaps(variosm, DMRs, type='within')
mcols(DMRs)<-DataFrame(aggregate(as.data.frame(mcols(variosm[idx@queryHits])), 
	by=list(DMR=idx@subjectHits), FUN=median))
DMRs


###################################################
### code chunk number 22: plotdmr
###################################################
s<-mcols(rhsmm@res)[,'SAMPLE']=='meth.diff' & seqnames(rhsmm@res) == 'chr1'
e<-seqnames(variosm) == 'chr1'
biomvRGviz(exprgr=variosm[e,'meth.diff'],  seggr=rhsmm@res[s], regionID='DMR', tofile=FALSE)


###################################################
### code chunk number 23: session
###################################################
sessionInfo()