### 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()