tutorial.R
In HelloRanges: Introduce *Ranges to bedtools users

### R code from vignette source 'tutorial.Rnw'

###################################################
### code chunk number 1: style-Sweave
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: loadData
###################################################
library(HelloRanges)
library(HelloRangesData)


###################################################
### code chunk number 3: setwd
###################################################
oldwd <- setwd(system.file("extdata", package="HelloRangesData"))


###################################################
### code chunk number 4: intersect-default
###################################################
code <- bedtools_intersect("-a cpg.bed -b exons.bed -g hg19.genome")
code


###################################################
### code chunk number 5: intersect-eval
###################################################
ans <- eval(code)
mcols(ans)$hit <- NULL
ans


###################################################
### code chunk number 6: intersect-simple
###################################################
code <- bedtools_intersect("-a cpg.bed -b exons.bed")
code


###################################################
### code chunk number 7: seqinfo
###################################################
genome <- eval(code[[2L]])
genome


###################################################
### code chunk number 8: intersect-genome
###################################################
bedtools_intersect("-a cpg.bed -b exons.bed -g hg19.genome")


###################################################
### code chunk number 9: intersect-genome-id
###################################################
bedtools_intersect("-a cpg.bed -b exons.bed -g hg19")


###################################################
### code chunk number 10: granges
###################################################
gr_a


###################################################
### code chunk number 11: rtracklayer-cpgs (eval = FALSE)
###################################################
## ucsc <- browserSession()
## genome(ucsc) <- "hg19"
## cpgs <- ucsc[["CpG Islands"]]


###################################################
### code chunk number 12: txdb-exons (eval = FALSE)
###################################################
## library(TxDb.Hsapiens.UCSC.hg19.knownGene)
## exons <- exons(TxDb.Hsapiens.UCSC.hg19.knownGene)


###################################################
### code chunk number 13: intersect-pairs
###################################################
pairs


###################################################
### code chunk number 14: export-pairs (eval = FALSE)
###################################################
## export(pairs, "pairs.bedpe")


###################################################
### code chunk number 15: intersect-ans
###################################################
ans


###################################################
### code chunk number 16: intersect-wa-wb
###################################################
bedtools_intersect("-a cpg.bed -b exons.bed -g hg19 -wa -wb")


###################################################
### code chunk number 17: intersect-wo
###################################################
bedtools_intersect("-a cpg.bed -b exons.bed -g hg19 -wo")


###################################################
### code chunk number 18: intersect-c
###################################################
bedtools_intersect("-a cpg.bed -b exons.bed -g hg19 -c")


###################################################
### code chunk number 19: intersect-v
###################################################
bedtools_intersect("-a cpg.bed -b exons.bed -g hg19 -v")


###################################################
### code chunk number 20: intersect-f
###################################################
bedtools_intersect("-a cpg.bed -b exons.bed -g hg19 -f 0.5 -wo")


###################################################
### code chunk number 21: intersect-performance
###################################################
a <- import("exons.bed")
b <- import("hesc.chromHmm.bed")
system.time(pintersect(findOverlapPairs(a, b, ignore.strand=TRUE),
                       ignore.strand=TRUE))


###################################################
### code chunk number 22: intersect-multiple
###################################################
code <- bedtools_intersect(
    paste("-a exons.bed",
          "-b cpg.bed,gwas.bed,hesc.chromHmm.bed -wa -wb -g hg19",
          "-names cpg,gwas,chromhmm"))
ans <- eval(code)
code


###################################################
### code chunk number 23: intersect-multiple-second
###################################################
second(ans)


###################################################
### code chunk number 24: merge
###################################################
bedtools_merge("-i exons.bed")


###################################################
### code chunk number 25: merge-count-overlaps
###################################################
code <- bedtools_merge("-i exons.bed -c 1 -o count")
code


###################################################
### code chunk number 26: merge-count-overlaps-ans
###################################################
eval(code)


###################################################
### code chunk number 27: merge-count-overlaps-direct
###################################################
identical(lengths(ans$grouping), ans$seqnames.count)


###################################################
### code chunk number 28: merge-close
###################################################
bedtools_merge("-i exons.bed -d 1000")


###################################################
### code chunk number 29: merge-multiple
###################################################
bedtools_merge("-i exons.bed -d 90 -c 1,4 -o count,collapse")


###################################################
### code chunk number 30: complement
###################################################
bedtools_complement("-i exons.bed -g hg19.genome")


###################################################
### code chunk number 31: genomecov-bg
###################################################
bedtools_genomecov("-i exons.bed -g hg19.genome -bg")


###################################################
### code chunk number 32: genomecov
###################################################
code <- bedtools_genomecov("-i exons.bed -g hg19.genome")
ans <- eval(code)
code


###################################################
### code chunk number 33: genomecov-cov
###################################################
cov


###################################################
### code chunk number 34: genomecov-ans
###################################################
ans


###################################################
### code chunk number 35: combine-genomecov
###################################################
code <- bedtools_genomecov("-i exons.bed -g hg19.genome -bga")
gr0 <- subset(eval(code), score == 0L) # compare to: awk '$4==0'
gr0


###################################################
### code chunk number 36: combine-intersect
###################################################
code <- R_bedtools_intersect("cpg.bed", gr0)
code


###################################################
### code chunk number 37: combine-pipe (eval = FALSE)
###################################################
## do_bedtools_genomecov("exons.bed", g="hg19.genome", bga=TRUE) %>% 
##     subset(score > 0L) %>%
##     do_bedtools_intersect("cpg.bed", .)    


###################################################
### code chunk number 38: coalescence-naive
###################################################
bedtools_coverage("-a cpg.bed -b exons.bed -hist -g hg19.genome")


###################################################
### code chunk number 39: coalescence-simplify
###################################################
genome <- import("hg19.genome")
gr_a <- import("cpg.bed", genome = genome)
gr_b <- import("exons.bed", genome = genome)
cov <- unname(coverage(gr_b)[gr_a])


###################################################
### code chunk number 40: coalescence-custom
###################################################
all_cov <- unlist(cov)
df <- as.data.frame(table(coverage=all_cov))
df$fraction <- df$Freq / length(all_cov)


###################################################
### code chunk number 41: tutorial.Rnw:615-616
###################################################
setwd(oldwd)


###################################################
### code chunk number 42: coalescence-plot
###################################################
plot((1-cumsum(fraction)) ~ as.integer(coverage), df, type="s",
     ylab = "fraction of bp > coverage", xlab="coverage")


###################################################
### code chunk number 43: tutorial.Rnw:622-623
###################################################
setwd(system.file("extdata", package="HelloRangesData"))


###################################################
### code chunk number 44: jaccard-example
###################################################
code <- bedtools_jaccard(
    paste("-a fHeart-DS16621.hotspot.twopass.fdr0.05.merge.bed",
          "-b fHeart-DS15839.hotspot.twopass.fdr0.05.merge.bed"))
heart_heart <- eval(code)
code <- bedtools_jaccard(
    paste("-a fHeart-DS16621.hotspot.twopass.fdr0.05.merge.bed",
          "-b fSkin_fibro_bicep_R-DS19745.hg19.hotspot.twopass.fdr0.05.merge.bed"))
heart_skin <- eval(code)
mstack(heart_heart=heart_heart, heart_skin=heart_skin)


###################################################
### code chunk number 45: jaccard-code
###################################################
code


###################################################
### code chunk number 46: jaccard-all
###################################################
files <- Sys.glob("*.merge.bed")
names(files) <- sub("\\..*", "", files)
ncores <- if (.Platform$OS.type == "windows") 1L else 4L
ans <- outer(files, files, 
             function(a, b) mcmapply(do_bedtools_jaccard, a, b, 
                                     mc.cores=ncores))
jaccard <- apply(ans, 1:2, function(x) x[[1]]$jaccard)


###################################################
### code chunk number 47: tutorial.Rnw:668-669
###################################################
setwd(oldwd)


###################################################
### code chunk number 48: jaccard-plot
###################################################
palette <- colorRampPalette(c("lightblue", "darkblue"))(9)
heatmap(jaccard, col=palette, margin=c(14, 14))


###################################################
### code chunk number 49: tutorial.Rnw:675-676
###################################################
setwd(system.file("extdata", package="HelloRangesData"))


###################################################
### code chunk number 50: answers-load
###################################################
genome <- import("hg19.genome")
exons <- import("exons.bed", genome=genome)
gwas <- import("gwas.bed", genome=genome)
hesc.chromHmm <- import("hesc.chromHmm.bed", genome=genome)


###################################################
### code chunk number 51: answer-1
###################################################
bedtools_complement("-i exons.bed -g hg19.genome")
## or without HelloRanges:
setdiff(as(seqinfo(exons), "GRanges"), unstrand(exons))


###################################################
### code chunk number 52: answer-2
###################################################
bedtools_closest("-a gwas.bed -b exons.bed -d")
## or 
distanceToNearest(gwas, exons)


###################################################
### code chunk number 53: answer-3
###################################################
code <- bedtools_makewindows("-g hg19.genome -w 500000")
code
windows <- unlist(eval(code))
R_bedtools_intersect(windows, exons, c=TRUE)
## or
str(countOverlaps(tileGenome(seqinfo(exons), tilewidth=500000), 
                  exons))


###################################################
### code chunk number 54: answer-4
###################################################
bedtools_intersect(
    paste("-a exons.bed -b <\"grep Enhancer hesc.chromHmm.bed\"",
          "-f 1.0 -wa -u"))
quote(length(ans))
## or
sum(exons %within% 
    subset(hesc.chromHmm, grepl("Enhancer", name)))


###################################################
### code chunk number 55: answer-5
###################################################
bedtools_intersect("-a gwas.bed -b exons.bed -u")
quote(length(gr_a)/length(ans))
## or
mean(gwas %over% exons)


###################################################
### code chunk number 56: answer-6
###################################################
bedtools_flank("-l 2 -r 2 -i exons.bed -g hg19.genome")
## or, bonus:
txid <- sub("_exon.*", "", exons$name)
tx <- split(exons, txid)
bounds <- range(tx)
transpliced <- lengths(bounds) > 1
introns <- unlist(psetdiff(unlist(bounds[!transpliced]), 
                           tx[!transpliced]))
Pairs(resize(introns, 2L), resize(introns, 2L, fix="end"))
## better way to get introns:
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
introns <- unlist(intronsByTranscript(txdb))


###################################################
### code chunk number 57: answer-7
###################################################
system.time(names(which.max(xtabs(width ~ name, 
                                  hesc.chromHmm))))
## or
names(which.max(sum(with(hesc.chromHmm, 
                         splitAsList(width, name)))))
## or
df <- aggregate(hesc.chromHmm, ~ name, totalWidth=sum(width))
df$name[which.max(df$totalWidth)]


###################################################
### code chunk number 58: restore-wd
###################################################
setwd(oldwd)