inst/doc/ctcfChipSeq.R

In igvR: igvR: integrative genomics viewer

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

## ---- eval=TRUE, echo=FALSE-------------------------------------------------------------------------------------------
knitr::include_graphics("igvR-ctcf-vignette-zoomedOut.png")

## ---- eval=TRUE, echo=FALSE-------------------------------------------------------------------------------------------
knitr::include_graphics("igvR-ctcf-vignette-zoomedIn.png")

## ----loadLibraries,  results='hide'-----------------------------------------------------------------------------------
#  library(igvR)
#  library(MotifDb)
#  library(AnnotationHub)
#  library(BSgenome.Hsapiens.UCSC.hg19)
#  library(phastCons100way.UCSC.hg19)
#  

## ----createLoad, results='hide'---------------------------------------------------------------------------------------
#  igv <- igvR()
#  setBrowserWindowTitle(igv, "CTCF ChIP-seq")
#  setGenome(igv, "hg19")

## ----initialDisplay,  results='hide'----------------------------------------------------------------------------------
#  showGenomicRegion(igv, "chr3:128,079,020-128,331,275")

## ----roi,  results='hide'---------------------------------------------------------------------------------------------
#  loc <- getGenomicRegion(igv)
#  which <- with(loc, GRanges(seqnames=chrom, ranges = IRanges(start, end)))
#  param <- ScanBamParam(which=which, what = scanBamWhat())
#  bamFile <- "~/github/ChIPseqMotifMatch/bulk/GSM749704/GSM749704_hg19_wgEncodeUwTfbsGm12878CtcfStdAlnRep1.bam"
#  file.exists(bamFile)
#  x <- readGAlignments(bamFile, use.names=TRUE, param=param)
#  track <- GenomicAlignmentTrack("ctcf bam", x, visibilityWindow=10000000, trackHeight=200)  # 30000 default
#  displayTrack(igv, track)

## ----narrow.peaks.track,  results='hide'------------------------------------------------------------------------------
#  filename <- system.file(package="igvR", "extdata", "GSM749704_hg19_chr19_peaks.narrowPeak")
#  tbl.pk <- get(load(system.file(package="igvR", "extdata", "GSM749704_hg19_chr19_subset.RData")))
#  track <- DataFrameQuantitativeTrack("macs2 peaks", tbl.pk, color="red", autoscale=TRUE)
#  displayTrack(igv, track)

## ----prepare.for.motif.match.and.display, results='hide'--------------------------------------------------------------
#  
#  dna <- with(loc, getSeq(BSgenome.Hsapiens.UCSC.hg19, chrom, start, end))
#  pfm.ctcf <- MotifDb::query(MotifDb, c("CTCF", "sapiens", "jaspar2018"), notStrings="ctcfl")
#  motif.name <- names(pfm.ctcf)[1]
#  pfm <- pfm.ctcf[[1]]

## ----motif.match, results='hide'--------------------------------------------------------------------------------------
#  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 <- loc$chrom
#  tbl.matches$start <- tbl.matches$start + loc$start
#  tbl.matches$end <- tbl.matches$end + loc$start

## ----prepare.motif.popup, results='hide'------------------------------------------------------------------------------
#  enableMotifLogoPopups(igv)
#  
#  tbl.matches$name <- paste0("MotifDb::", motif.name)
#  tbl.matches <- tbl.matches[, c("chrom", "start", "end", "name", "score")]
#  dim(tbl.matches)

## ----display.both.annotation.and.quantitative.tracks, results='hide'--------------------------------------------------
#  track <- DataFrameAnnotationTrack("CTCF-MA0139.1", tbl.matches[, 1:4], color="random",
#                                    trackHeight=25)
#  displayTrack(igv, track)
#  
#  track <- DataFrameQuantitativeTrack("MA0139.1 scores", tbl.matches[, c(1,2,3,5)],
#                                       color="random", autoscale=FALSE, min=8, max=12)
#  displayTrack(igv, track)
#  

## ----conservation, results='hide'-------------------------------------------------------------------------------------
#  loc <- getGenomicRegion(igv)
#  starts <- with(loc, seq(start, end, by=5))
#  ends <- starts + 5
#  count <- length(starts)
#  tbl.blocks <- data.frame(chrom=rep(loc$chrom, count), start=starts, end=ends, stringsAsFactors=FALSE)
#  dim(tbl.blocks)
#  
#  tbl.cons100 <- as.data.frame(gscores(phastCons100way.UCSC.hg19, GRanges(tbl.blocks)), stringsAsFactors=FALSE)
#  tbl.cons100$chrom <- as.character(tbl.cons100$seqnames)
#  tbl.cons100 <- tbl.cons100[, c("chrom", "start", "end", "default")]
#  track <- DataFrameQuantitativeTrack("cons100", tbl.cons100, autoscale=TRUE, color="darkgreen")
#  displayTrack(igv, track)

## ----methylationTracks,  results='hide'-------------------------------------------------------------------------------
#  ah <- AnnotationHub()
#  ah.human <- subset(ah, species == "Homo sapiens")
#  
#  histone.tracks <- AnnotationHub::query(ah.human, c("H3K4me3", "Gm12878", "Peak", "narrow"))  # 3 tracks
#  descriptions <- histone.tracks$description
#  titles <- histone.tracks$title
#  colors <- rep(terrain.colors(6), 4)
#  color.index <- 0
#  
#  roi <- with(getGenomicRegion(igv), GRanges(seqnames=chrom, IRanges(start=start, end=end)))
#  
#  for(i in seq_len(length(histone.tracks))){
#     name <- names(histone.tracks)[i]
#     color.index <- color.index + 1
#     gr <- histone.tracks[[name]]
#     ov <- findOverlaps(gr, roi)
#     histones <- gr[queryHits(ov)]
#     track.histones <- GRangesQuantitativeTrack(titles[i], histones[, "pValue"],
#                                                color=colors[color.index], trackHeight=50,
#                                                autoscale=TRUE)
#     displayTrack(igv, track.histones)
#     Sys.sleep(5)
#     } # for track
#  

## ----queryAHforPromoters,  results='hide'-----------------------------------------------------------------------------
#  
#  

## ----zoom.in, results='hide'------------------------------------------------------------------------------------------
#  showGenomicRegion(igv, "chr3:128,216,201-128,217,324")

## ----sessionInfo------------------------------------------------------------------------------------------------------
#  sessionInfo()