## ---- include=FALSE----------------------------------------------------------- library(knitr) library(regioneR) opts_chunk$set(concordance=FALSE) set.seed(21666641) ## ----------------------------------------------------------------------------- A <- createRandomRegions(nregions=50, length.mean=5000000, length.sd=3000000) B <- c(A[1:25], createRandomRegions(nregions=25, length.mean=500000, length.sd=30000)) numOverlaps(A, B, count.once=TRUE) numOverlaps(randomizeRegions(A), B, count.once=TRUE) ## ----------------------------------------------------------------------------- pt <- overlapPermTest(A=A, B=B, ntimes=50) pt ## ---- fig.height=6,fig.width=8------------------------------------------------ plot(pt) ## ---- fig.height=6,fig.width=8------------------------------------------------ lz <- localZScore(pt=pt, A=A, B=B) plot(lz) ## ----------------------------------------------------------------------------- special <- toGRanges(system.file("extdata", "my.special.genes.txt", package="regioneR")) all.genes <- toGRanges(system.file("extdata", "all.genes.txt", package="regioneR")) altered <- toGRanges(system.file("extdata", "my.altered.regions.txt", package="regioneR")) length(special) length(all.genes) length(altered) ## ----------------------------------------------------------------------------- numOverlaps(special, altered) ## ----------------------------------------------------------------------------- random.RS <- resampleRegions(special, universe=all.genes) random.RS ## ----------------------------------------------------------------------------- random.RS <- resampleRegions(special, universe=all.genes) numOverlaps(random.RS, altered) random.RS <- resampleRegions(special, universe=all.genes) numOverlaps(random.RS, altered) random.RS <- resampleRegions(special, universe=all.genes) numOverlaps(random.RS, altered) random.RS <- resampleRegions(special, universe=all.genes) numOverlaps(random.RS, altered) ## ---- eval=FALSE-------------------------------------------------------------- # # NOT RUN # pt <- permTest(A=my.regions, B=repeats, randomize.function=randomizeRegions, # evaluate.function=numOverlaps) ## ---- eval=FALSE-------------------------------------------------------------- # # NOT RUN # pt <- permTest(A=my.genes, randomize.function=resampleRegions, universe=all.genes, # evaluate.function=meanInRegions, x=methylation.levels.450K) ## ----------------------------------------------------------------------------- pt <- permTest(A=special, ntimes=50, randomize.function=resampleRegions, universe=all.genes, evaluate.function=numOverlaps, B=altered, verbose=FALSE) ## ----------------------------------------------------------------------------- pt summary(pt) ## ---- fig.height=6,fig.width=8------------------------------------------------ plot(pt) ## ---- fig.height=6,fig.width=8------------------------------------------------ regular <- toGRanges(system.file("extdata", "my.regular.genes.txt", package="regioneR")) length(regular) numOverlaps(regular, altered) pt.reg <- permTest(A=regular, ntimes=50, randomize.function=resampleRegions, universe=all.genes, evaluate.function=numOverlaps, B=altered, verbose=FALSE) pt.reg plot(pt.reg) ## ---- eval=FALSE-------------------------------------------------------------- # #NOT RUN - See Figure 1 # pt.5000 <- permTest(A=special, ntimes=5000, randomize.function=resampleRegions, # universe=all.genes, evaluate.function=numOverlaps, B=altered, verbose=TRUE) # plot(pt.5000) ## ---- eval=FALSE-------------------------------------------------------------- # #NOT RUN - See Figure 2 # pt.5000.reg <- permTest(A=regular, ntimes=5000, randomize.function=resampleRegions, # universe=all.genes, evaluate.function=numOverlaps, B=altered, verbose=TRUE) # plot(pt.5000.reg) ## ----------------------------------------------------------------------------- A <- toGRanges(data.frame(chr=c("chr1", "chr1", "chr1"), start=c(20000, 50000, 100000), end=c(22000, 70000, 400000))) ## ----------------------------------------------------------------------------- randomizeRegions(A, genome="hg19") ## ----------------------------------------------------------------------------- randomizeRegions(A, genome="hg19", per.chromosome=TRUE) ## ----------------------------------------------------------------------------- circularRandomizeRegions(A, genome="hg19", mask=NA) ## ----------------------------------------------------------------------------- gcContent <- function(A, bsgenome, ...) { A <- toGRanges(A) reg.seqs <- getSeq(bsgenome, A) base.frequency <- alphabetFrequency(reg.seqs) gc.pct <- (sum(base.frequency[,"C"]) + sum(base.frequency[,"G"]))/sum(width(A)) return(gc.pct) } ## ----------------------------------------------------------------------------- permuteRegionsMetadata <- function(A, ...) { A <- toGRanges(A) mcols(A)[,1] <- mcols(A)[sample(length(A)),1] return(A) } ## ----------------------------------------------------------------------------- A <- createRandomRegions(nregions=50, length.mean=5000000, length.sd=3000000) B <- c(A[1:25], createRandomRegions(nregions=25, length.mean=500000, length.sd=30000)) #Without mc.set.seed=FALSE set.seed(123) overlapPermTest(A, B, ntimes=10, force.parallel=TRUE) set.seed(123) overlapPermTest(A, B, ntimes=10, force.parallel=TRUE) #With mc.set.seed=FALSE set.seed(123) overlapPermTest(A, B, ntimes=10, mc.set.seed=FALSE, force.parallel=TRUE) set.seed(123) overlapPermTest(A, B, ntimes=10, mc.set.seed=FALSE, force.parallel=TRUE) ## ---- fig.height=6,fig.width=8------------------------------------------------ pt <- permTest(A=special, ntimes=50, randomize.function=resampleRegions, universe=all.genes, evaluate.function=numOverlaps, B=altered, verbose=FALSE) lz <- localZScore(A=special, pt=pt, window=10*mean(width(special)), step=mean(width(special))/2, B=altered) plot(lz) ## ---- fig.height=6,fig.width=8------------------------------------------------ genome <- filterChromosomes(getGenome("hg19"), keep.chr="chr1") B <- createRandomRegions(nregions=100, length.mean=10000, length.sd=20000, genome=genome, non.overlapping=FALSE) A <- B[sample(20)] pt <- overlapPermTest(A=A, B=B, ntimes=100, genome=genome, non.overlapping=FALSE) pt lz <- localZScore(A=A, B=B, pt=pt, window=10*mean(width(A)), step=mean(width(A))/2 ) plot(lz) ## ----------------------------------------------------------------------------- A <- data.frame(chr=1, start=c(1,15,24), end=c(10,20,30), x=c(1,2,3), y=c("a","b","c")) B <- toGRanges(A) B ## ----------------------------------------------------------------------------- toDataframe(B) ## ----------------------------------------------------------------------------- human.genome <- getGenomeAndMask("hg19", mask=NA)$genome human.canonical <- filterChromosomes(human.genome, organism="hg") listChrTypes() human.autosomal <- filterChromosomes(human.genome, organism="hg", chr.type="autosomal") human.123 <- filterChromosomes(human.genome, keep.chr=c("chr1", "chr2", "chr3")) ## ----------------------------------------------------------------------------- overlaps.df <- overlapRegions(A, B) overlaps.df ## ----------------------------------------------------------------------------- overlaps.df <- overlapRegions(A, B, type="BinA", get.pctA=TRUE) overlaps.df ## ----------------------------------------------------------------------------- overlaps.df <- overlapRegions(A, B, min.bases=5) overlaps.df ## ----------------------------------------------------------------------------- overlaps.bool <- overlapRegions(A, B, min.bases=5, only.boolean=TRUE) overlaps.bool ## ----------------------------------------------------------------------------- overlaps.int <- overlapRegions(A, B, min.bases=5, only.count=TRUE) overlaps.int ## ---- eval=FALSE-------------------------------------------------------------- # #NOT RUN # set.seed(12345) # cpgHMM <- toGRanges("http://www.haowulab.org/software/makeCGI/model-based-cpg-islands-hg19.txt") # promoters <- toGRanges("http://gattaca.imppc.org/regioner/data/UCSC.promoters.hg19.bed") ## ---- eval=FALSE-------------------------------------------------------------- # #NOT RUN # cpgHMM <- filterChromosomes(cpgHMM, organism="hg", chr.type="canonical") # promoters <- filterChromosomes(promoters, organism="hg", chr.type="canonical") ## ---- eval=FALSE-------------------------------------------------------------- # #NOT RUN # cpgHMM_2K <- sample(cpgHMM, 2000) # # pt <- overlapPermTest(cpgHMM_2K, promoters, ntimes=1000, genome="hg19", count.once=TRUE) # pt # P-value: 0.000999000999000999 # Z-score: 60.5237 # Number of iterations: 1000 # Alternative: greater # Evaluation of the original region set: 614 # Evaluation function: numOverlaps # Randomization function: randomizeRegions # mean(pt$permuted) # 79.087 ## ---- eval=FALSE-------------------------------------------------------------- # #NOT RUN # plot(pt) ## ---- eval=FALSE-------------------------------------------------------------- # #NOT RUN # set.seed(12345) # download.file("http://hgdownload-test.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsSydhHepg2Rad21IggrabUniPk.narrowPeak.gz", "Rad21.gz") # # download.file("http://hgdownload-test.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsUwHepg2CtcfUniPk.narrowPeak.gz", "Ctcf.gz") # # HepG2_Rad21 <- toGRanges(gzfile("Rad21.gz"), header=FALSE) # HepG2_Ctcf <- toGRanges(gzfile("Ctcf.gz"), header=FALSE) # # promoters <- toGRanges("http://gattaca.imppc.org/regioner/data/UCSC.promoters.hg19.bed") ## ---- eval=FALSE-------------------------------------------------------------- # #NOT RUN # promoters <- filterChromosomes(promoters, organism="hg19") # HepG2_Rad21_5K <- sample(HepG2_Rad21, 5000) ## ---- eval=FALSE-------------------------------------------------------------- # #NOT RUN # pt_Rad21_5k_vs_Ctcf <- permTest(A=HepG2_Rad21_5K, B=HepG2_Ctcf, ntimes=1000, # randomize.function=circularRandomizeRegions, # evaluate.function=numOverlaps, count.once=TRUE, # genome="hg19", mc.set.seed=FALSE, mc.cores=4) # # pt_Rad21_5k_vs_Prom <- permTest(A=HepG2_Rad21_5K, B=promoters, ntimes=1000, # randomize.function=circularRandomizeRegions, # evaluate.function=numOverlaps, count.once=TRUE, # genome="hg19", mc.set.seed=FALSE, mc.cores=4) # # pt_Rad21_5k_vs_Ctcf # # pt_Rad21_5k_vs_Prom # # plot(pt_Rad21_5k_vs_Ctcf, main="Rad21_5K vs CTCF") # plot(pt_Rad21_5k_vs_Prom, main="Rad21_5K vs Promoters") ## ---- eval=FALSE-------------------------------------------------------------- # #NOT RUN # lz_Rad21_vs_Ctcf_1 <- localZScore(A=HepG2_Rad21_5K, B=HepG2_Ctcf, # pt=pt_Rad21_5k_vs_Ctcf, # window=1000, step=50, count.once=TRUE) # # lz_Rad21_vs_Prom_1 <- localZScore(A=HepG2_Rad21_5K, B=promoters, # pt=pt_Rad21_5k_vs_Prom, # window=1000, step=50, count.once=TRUE) # # plot(lz_Rad21_vs_Ctcf_1, main="Rad21_5k vs CTCF (1Kbp)") # plot(lz_Rad21_vs_Prom_1, main="Rad21_5k vs promoters (1Kbp)") # ## ---- eval=FALSE-------------------------------------------------------------- # #NOT RUN # lz_Rad21_vs_Ctcf_2 <- localZScore(A=HepG2_Rad21_5K, B=HepG2_Ctcf, # pt=pt_Rad21_5k_vs_Ctcf, # window=10000, step=500, count.once=TRUE) # # lz_Rad21_vs_Prom_2 <- localZScore(A=HepG2_Rad21_5K, B=promoters, # pt=pt_Rad21_5k_vs_Prom, # window=10000, step=500, count.once=TRUE) # # plot(lz_Rad21_vs_Ctcf_2, main="Rad21_5k vs CTCF (10Kbp)") # plot(lz_Rad21_vs_Prom_2, main="Rad21_5k vs promoters (10Kbp)") ## ---- sessionInfo------------------------------------------------------------- sessionInfo()