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

###################################################
### chunk number 1: setup
###################################################
options(width=50)


###################################################
### chunk number 2: chipseq-cstest
###################################################
library(chipseq)
data(cstest)
names(cstest)
head(cstest$ctcf, 1)


###################################################
### chunk number 3: chipseq-estimate.mean.fraglen eval=FALSE
###################################################
## fraglen <- 
##   estimate.mean.fraglen(cstest$ctcf, 
##                         method = "correlation",
##                         seqLen = 35)


###################################################
### chunk number 4: chipseq-estimate.mean.fraglen
###################################################
library(lazycache)
cache(fraglen = 
      estimate.mean.fraglen(cstest$ctcf, 
                            method = "correlation", 
                            seqLen = 35))


###################################################
### chunk number 5: chipseq-fraglen
###################################################
fraglen
median(fraglen)


###################################################
### chunk number 6: chipseq-resize
###################################################
ctcf.ext <- resize(cstest$ctcf, 
                   width = median(fraglen))
head(ctcf.ext, 3)


###################################################
### chunk number 7: chipseq-coverage
###################################################
cov.ctcf <- coverage(ctcf.ext)
cov.ctcf$chr10


###################################################
### chunk number 8: chipseq-slice
###################################################
islands <- slice(cov.ctcf, lower = 1)
head(islands$chr10, 3)


###################################################
### chunk number 9: chipseq-slice-peaks
###################################################
peak_viewsList <- slice(cov.ctcf, lower = 8)
peak_rangesList <- ranges(peak_viewsList)
peaks <- as(peak_rangesList, "GRanges")
head(peaks, 3)


###################################################
### chunk number 10: findPeaks
###################################################
findPeaks <- function(reads) {
  fraglen <- estimate.mean.fraglen(reads, 
                                   method = "correlation", 
                                   seqLen = 35)
  reads_ext <- resize(reads, fraglen)
  cov <- coverage(reads_ext)
  slice(cov, lower = 8)
}


###################################################
### chunk number 11: chipseq-viewMaxs
###################################################
values(peaks)$max <- 
  unlist(viewMaxs(peak_viewsList))


###################################################
### chunk number 12: chipseq-viewSums
###################################################
values(peaks)$sum <- 
  unlist(viewSums(peak_viewsList))


###################################################
### chunk number 13: chipseq-viewRangeMaxs
###################################################
values(peaks)$summits <- 
  unlist(viewRangeMaxs(peak_viewsList))


###################################################
### chunk number 14: chipseq-promoters
###################################################
library(TxDb.Mmusculus.UCSC.mm9.knownGene)
mm9_tx <- transcripts(Mmusculus_UCSC_mm9_knownGene_TxDb,
                      columns = "gene_id")


###################################################
### chunk number 15: chipseq-gene-id
###################################################
gene_id <- values(mm9_tx)$gene_id
all(elementLengths(gene_id) <= 1)
flat_gene_id <- character(length(mm9_tx))
flat_gene_id[elementLengths(gene_id) == 1] <- 
  unlist(gene_id)
values(mm9_tx)$gene_id <- flat_gene_id


###################################################
### chunk number 16: chipseq-promoters
###################################################
promoters <- resize(flank(mm9_tx, 500), 700)


###################################################
### chunk number 17: chipseq-promoters-in
###################################################
values(peaks)$in_promoter <- peaks %in% promoters
table(values(peaks)$in_promoter)


###################################################
### chunk number 18: chipseq-promoters-match
###################################################
values(peaks)$in_promoter_of <- 
  values(promoters)$gene_id[match(peaks, promoters)]
head(subset(values(peaks), in_promoter)$in_promoter_of, 3)


###################################################
### chunk number 19: chipseq-all-peaks eval=FALSE
###################################################
## ctcf_peaks <- findPeaks(cstest$ctcf)
## gfp_peaks <- findPeaks(cstest$gfp)


###################################################
### chunk number 20: chipseq-all-peaks-real
###################################################
lazycache::cache(ctcf_peaks = findPeaks(cstest$ctcf))
lazycache::cache(gfp_peaks = findPeaks(cstest$gfp))


###################################################
### chunk number 21: chipseq-diffPeakSummary
###################################################
peak_summary <- diffPeakSummary(ctcf_peaks, gfp_peaks)
colnames(peak_summary)


###################################################
### chunk number 22: chipseq-SummarizedExperiment
###################################################
peak_summary_gr <- as(ranges(peak_summary), "GRanges")
max_matrix <- with(peak_summary, cbind(maxs1, maxs2))
SummarizedExperiment(max_matrix, rowData = peak_summary_gr)


###################################################
### chunk number 23: rnaseq-readGappedAlignments
###################################################
library(leeBamViews)
bams <- dir(system.file("bam", package="leeBamViews"), 
            full = TRUE, pattern = "bam$")
reads_ga <- readGappedAlignments(bams[1])
head(reads_ga, 1)


###################################################
### chunk number 24: rnaseq-grglist
###################################################
reads_grl <- grglist(reads_ga)


###################################################
### chunk number 25: rnaseq-transcripts
###################################################
data(leeUnn)
leeUnn <- 
  subset(leeUnn, lengthWithoutMask > 0 & !is.na(chr))
leeUnn$strand <- c("-", "*", "+")[leeUnn$strand + 2]
sc2_tx <- with(leeUnn, 
               GRanges(sprintf("Scchr%02d", chr), 
                       IRanges(start, end), strand))
seqlevels(sc2_tx)[length(seqlevels(sc2_tx))] <- "Scmito"


###################################################
### chunk number 26: rnaseq-count-exons
###################################################
values(sc2_tx)$counts <- 
  countOverlaps(sc2_tx, reads_grl, ignore.strand = TRUE)


###################################################
### chunk number 27: rnaseq-junctions
###################################################
ol <- findOverlaps(reads_grl, sc2_tx, 
                   ignore.strand = TRUE)
reads_factor <- factor(queryHits(ol), 
                       seq_len(length(reads_grl)))
values(reads_grl)$tx_hits <- 
  split(subjectHits(ol), reads_factor)
head(table(elementLengths(values(reads_grl)$exon_hits)))