## ----style-knitr, eval=TRUE, echo=FALSE, results="asis"------------------
BiocStyle::latex()


## ----include=FALSE-------------------------------------------------------
library(knitr)
library(regioneR)
opts_chunk$set(concordance=FALSE)


## ------------------------------------------------------------------------
special <- toGRanges("http://gattaca.imppc.org/regioner/data/my.special.genes.bed")
all.genes <- toGRanges("http://gattaca.imppc.org/regioner/data/all.genes.bed")
altered <- toGRanges("http://gattaca.imppc.org/regioner/data/my.altered.regions.bed")
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=overlapRegions)


## ----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("http://gattaca.imppc.org/regioner/data/my.regular.genes.bed")
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")


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


## ----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")$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