### R code from vignette source 'vignettes/GenomicRanges/inst/doc/countFeatureHits.Rnw'

###################################################
### code chunk number 1: firstExample (eval = FALSE)
###################################################
## library(TxDb.Hsapiens.UCSC.hg18.knownGene)
## txdb <- Hsapiens_UCSC_hg18_knownGene_TxDb
## ranges <- exonsBy(txdb, "tx")
## bamFiles <- c("sample1.bam", "sample2.bam")
## 
## countsTable <-
##     as.data.frame(tx_id=rep(names(ranges), elementLengths(ranges)),
##                   exon_id=values(unlist(ranges))[["exon_id"]],
##                   assays(countFeatureHits(bamFiles, ranges))$counts
##                   )


###################################################
### code chunk number 2: simple
###################################################
library(GenomicRanges)
rd <- GappedAlignments("a", rname = Rle("chr1"), pos = as.integer(100),
    cigar = "300M", strand = strand("+"))

gr1 <- GRanges("chr1", IRanges(start=50, width=150), strand="+")
gr2 <- GRanges("chr1", IRanges(start=350, width=150), strand="+")
gr <- c(gr1, gr2)
grl <- GRangesList(c(gr1, gr2))


###################################################
### code chunk number 3: simpleGRanges
###################################################
data.frame(union = assays(countFeatureHits(rd, gr))$counts,
           intStrict = assays(countFeatureHits(rd, gr,
               mode="IntersectionStrict"))$counts,
           intNotEmpty = assays(countFeatureHits(rd, gr,
               mode="IntersectionNotEmpty"))$counts)


###################################################
### code chunk number 4: simpleGRangesList
###################################################
data.frame(union = assays(countFeatureHits(rd, grl))$counts,
           intStrict = assays(countFeatureHits(rd, grl,
               mode="IntersectionStrict"))$counts,
           intNotEmpty = assays(countFeatureHits(rd, grl,
               mode="IntersectionNotEmpty"))$counts)


###################################################
### code chunk number 5: data
###################################################
group_id <- c("A", "B", "C", "C", "D", "D", "E", "F", "G", "G", "H", "H")
features <- GRanges(
    seqnames = Rle(c("chr1", "chr2", "chr1", "chr1", "chr2", "chr2",
        "chr1", "chr1", "chr2", "chr2", "chr1", "chr1")),
    strand = strand(rep("+", length(group_id))),
    ranges = IRanges(
        start=c(1000, 2000, 3000, 3600, 7000, 7500, 4000, 4000, 3000, 3350, 5000, 5400),
        width=c(500, 900, 500, 300, 600, 300, 500, 900, 150, 200, 500, 500)),
   DataFrame(group_id)
)

reads <- GappedAlignments(
    names = c("a","b","c","d","e","f","g"),
    rname = Rle(c(rep(c("chr1", "chr2"), 3), "chr1")),
    pos = as.integer(c(1400, 2700, 3400, 7100, 4000, 3100, 5200)),
    cigar = c("500M", "100M", "300M", "500M", "300M", "50M200N50M", "50M150N50M"),
    strand = strand(rep.int("+", 7L)))



###################################################
### code chunk number 6: GRanges
###################################################
data.frame(union = assays(countFeatureHits(reads, features))$counts,
           intStrict = assays(countFeatureHits(reads, features,
               mode="IntersectionStrict"))$counts,
           intNotEmpty = assays(countFeatureHits(reads, features,
               mode="IntersectionNotEmpty"))$counts)


###################################################
### code chunk number 7: lst
###################################################
lst <- split(features, values(features)[["group_id"]])
length(lst)



###################################################
### code chunk number 8: GRangesList
###################################################
data.frame(union = assays(countFeatureHits(reads, lst))$counts,
           intStrict = assays(countFeatureHits(reads, lst,
               mode="IntersectionStrict"))$counts,
           intNotEmpty = assays(countFeatureHits(reads, lst,
               mode="IntersectionNotEmpty"))$counts)