## ----setup, include = FALSE-------------------------------------------------------------------------------------------
options(width=120)
knitr::opts_chunk$set(
   collapse = TRUE,
   eval=interactive(),
   echo=TRUE,
   comment = "#>"
)

## ----eval=FALSE-------------------------------------------------------------------------------------------------------
#  library(igvR)
#  igv <- igvR()
#  setBrowserWindowTitle(igv, "CTCF ChIP-seq")
#  setGenome(igv, "hg19")
#  showGenomicRegion(igv, "chr3:128,079,020-128,331,275")
#     # or
#  showGenomicRegion(igv, "GATA2")
#  for(i in 1:4) zoomOut(igv)

## ---- eval=TRUE, echo=FALSE-------------------------------------------------------------------------------------------
knitr::include_graphics("images/ctcfBam-01.png")

## ----eval=FALSE-------------------------------------------------------------------------------------------------------
#  roi <- getGenomicRegion(igv)
#  gr.roi <- with(roi, GRanges(seqnames=chrom, ranges = IRanges(start, end)))
#  param <- ScanBamParam(which=gr.roi, what = scanBamWhat())
#  bamFile <- system.file(package="igvR", "extdata", "ctcf-gata2", "gata2-region-hg19.bam")
#  alignments <- readGAlignments(bamFile, use.names=TRUE, param=param)

## ----eval=FALSE-------------------------------------------------------------------------------------------------------
#  track <- GenomicAlignmentTrack(trackName="ctcf bam", alignments, visibilityWindow=10000000, trackHeight=200)
#  displayTrack(igv, track)

## ---- eval=TRUE, echo=FALSE-------------------------------------------------------------------------------------------
knitr::include_graphics("images/ctcfBam-02.png")

## ----eval=FALSE-------------------------------------------------------------------------------------------------------
#  lapply(tbl.pk, class)
#  head(tbl.pk)
#  
#  constructed from a data.frame.
#  
#  

## ----eval=FALSE-------------------------------------------------------------------------------------------------------
#  narrowPeaksFile <- system.file(package="igvR", "extdata", "ctcf-gata2",
#                                 "gata2-region-macs2-narrowPeaks.RData")
#  tbl.pk <- get(load(narrowPeaksFile))
#  dim(tbl.pk) # 109 4
#  head(tbl.pk)
#    #      chrom     start       end score
#    # 6381 chr10 127910682 127910864    27
#    # 6382 chr10 128075644 128075811    89
#    # 6383 chr10 128259852 128259984    27
#    # 6384 chr10 128286655 128286920    78
#    # 6385 chr10 128437706 128437938    89
#    # 8827 chr11 127965327 127965489    70
#  
#  unlist(lapply(tbl.pk, class))
#  #        chrom       start         end       score
#  #  "character"   "integer"   "integer"   "integer"
#  
#  track <- DataFrameQuantitativeTrack("macs2 peaks", tbl.pk, color="red", autoscale=TRUE)
#  displayTrack(igv, track)
#  

## ---- eval=TRUE, echo=FALSE-------------------------------------------------------------------------------------------
knitr::include_graphics("images/ctcfBam-03.png")

## ---- out.width="500px", eval=TRUE, echo=FALSE------------------------------------------------------------------------
knitr::include_graphics("images/ctcfMotif.png")

## ----eval=FALSE-------------------------------------------------------------------------------------------------------
#     # get the DNA sequence in the current region
#  dna <- with(roi, getSeq(BSgenome.Hsapiens.UCSC.hg19, chrom, start, end))
#     # get the first of three motifs for CTCF.  (a more thorough study would use all three)
#  pfm.ctcf <- MotifDb::query(MotifDb, c("CTCF", "sapiens", "jaspar2022"), notStrings="ctcfl")
#  motif.name <- names(pfm.ctcf)[1]
#  pfm <- pfm.ctcf[[1]]
#  
#     # Find matches >= 80% of this motif in the sequence.  create a suitable
#     # data.frame for another DataFrameQuantitativeTrack
#  
#  hits.forward <- matchPWM(pfm, as.character(dna), with.score=TRUE, min.score="80%")
#  hits.reverse <- matchPWM(reverseComplement(pfm), as.character(dna), with.score=TRUE, min.score="80%")
#  
#  tbl.forward <- as.data.frame(ranges(hits.forward))
#  tbl.reverse <- as.data.frame(ranges(hits.reverse))
#  tbl.forward$score <- mcols(hits.forward)$score
#  tbl.reverse$score <- mcols(hits.reverse)$score
#  
#  tbl.matches <- rbind(tbl.forward, tbl.reverse)
#  tbl.matches$chrom <- roi$chrom
#  tbl.matches$start <- tbl.matches$start + roi$start
#  
#  tbl.matches$end <- tbl.matches$end + roi$start
#  
#  tbl.matches$name <- paste0("MotifDb::", motif.name)
#  tbl.matches <- tbl.matches[, c("chrom", "start", "end", "name", "score")]
#  dim(tbl.matches) # 25 5
#  head(tbl.matches)
#      #   chrom     start       end                                       name    score
#      # 1  chr3 128110910 128110928 MotifDb::Hsapiens-jaspar2018-CTCF-MA0139.1 10.70369
#      # 2  chr3 128114573 128114591 MotifDb::Hsapiens-jaspar2018-CTCF-MA0139.1 10.36891
#      # 3  chr3 128127658 128127676 MotifDb::Hsapiens-jaspar2018-CTCF-MA0139.1 10.44666
#      # 4  chr3 128138376 128138394 MotifDb::Hsapiens-jaspar2018-CTCF-MA0139.1 10.54861
#      # 5  chr3 128139280 128139298 MotifDb::Hsapiens-jaspar2018-CTCF-MA0139.1 10.51689
#      # 6  chr3 128173128 128173146 MotifDb::Hsapiens-jaspar2018-CTCF-MA0139.1 10.37987
#  

## ----eval=FALSE-------------------------------------------------------------------------------------------------------
#  track <- DataFrameQuantitativeTrack("MA0139.1 scores", tbl.matches[, c(1,2,3,5)],
#                                      color="random", autoscale=FALSE, min=8, max=12)
#  displayTrack(igv, track)

## ---- eval=TRUE, echo=FALSE-------------------------------------------------------------------------------------------
knitr::include_graphics("images/ctcfBam-04.png")

## ----eval=FALSE-------------------------------------------------------------------------------------------------------
#  
#  roi <- getGenomicRegion(igv)
#  bigwig.file <- system.file(package="igvR", "extdata", "ctcf-gata2", "gata2-region-h3k4me3.bw")
#  
#  bw.slice <- import(bigwig.file, which=gr.roi)
#  track <- GRangesQuantitativeTrack("h3k4me3", bw.slice, autoscale=TRUE)
#  displayTrack(igv, track)
#  

## ---- eval=TRUE, echo=FALSE-------------------------------------------------------------------------------------------
knitr::include_graphics("images/ctcfBam-05.png")

## ----eval=FALSE-------------------------------------------------------------------------------------------------------
#  showGenomicRegion(igv, "chr3:128,202,505-128,222,868")

## ---- eval=TRUE, echo=FALSE-------------------------------------------------------------------------------------------
knitr::include_graphics("images/ctcfBam-07.png")

## ----eval=TRUE--------------------------------------------------------------------------------------------------------
sessionInfo()