## ----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()