## ---- echo=FALSE, results="hide", warning=FALSE------------------------------- suppressPackageStartupMessages({ library(ChIPpeakAnno) library(EnsDb.Hsapiens.v75) library(TxDb.Hsapiens.UCSC.hg19.knownGene) library(org.Hs.eg.db) }) knitr::opts_chunk$set(warning=FALSE, message=FALSE) ## ----quickStart--------------------------------------------------------------- ## First, load the ChIPpeakAnno package library(ChIPpeakAnno) ## ----import------------------------------------------------------------------- path <- system.file("extdata", "Tead4.broadPeak", package="ChIPpeakAnno") peaks <- toGRanges(path, format="broadPeak") peaks[1:2] ## ----annotationData----------------------------------------------------------- library(EnsDb.Hsapiens.v75) annoData <- toGRanges(EnsDb.Hsapiens.v75) annoData[1:2] ## ----annotate----------------------------------------------------------------- ## keep the seqnames in the same style seqlevelsStyle(peaks) <- seqlevelsStyle(annoData) ## do annotation by nearest TSS anno <- annotatePeakInBatch(peaks, AnnotationData=annoData) anno[1:2] # A pie chart can be used to demonstrate the overlap features of the peaks. pie1(table(anno$insideFeature)) ## ----addIDs------------------------------------------------------------------- library(org.Hs.eg.db) anno <- addGeneIDs(anno, orgAnn="org.Hs.eg.db", feature_id_type="ensembl_gene_id", IDs2Add=c("symbol")) head(anno) ## ----------------------------------------------------------------------------- library(TxDb.Hsapiens.UCSC.hg19.knownGene) annoData <- toGRanges(TxDb.Hsapiens.UCSC.hg19.knownGene) annoData[1:2] seqlevelsStyle(peaks) <- seqlevelsStyle(annoData) ## ----------------------------------------------------------------------------- anno <- annotatePeakInBatch(peaks, AnnotationData=annoData, output="overlapping", FeatureLocForDistance="TSS", bindingRegion=c(-2000, 300)) anno$symbol <- xget(anno$feature, org.Hs.egSYMBOL) head(anno) ## ---- fig.height=3, fig.width=8----------------------------------------------- anno <- annotatePeakInBatch(peaks, AnnotationData=annoData, output="nearestBiDirectionalPromoters", bindingRegion=c(-5000, 500)) anno$symbol <- xget(anno$feature, org.Hs.egSYMBOL) anno[anno$peak=="peak12725"] ## ----trackViewer-------------------------------------------------------------- library(trackViewer) gr <- peak <- peaks["peak12725"] start(gr) <- start(gr) - 5000 end(gr) <- end(gr) + 5000 if(.Platform$OS.type != "windows"){ peak12725 <- importScore(file=system.file("extdata", "Tead4.bigWig", package="ChIPpeakAnno"), ranges=peak, format = "BigWig") }else{## rtracklayer can not import bigWig files on Windows load(file.path(dirname(path), "cvglist.rds")) peak12725 <- Views(cvglists[["Tead4"]][[as.character(seqnames(peak))]], start(peak), end(peak)) peak12725 <- viewApply(peak12725, as.numeric) tmp <- rep(peak, width(peak)) width(tmp) <- 1 tmp <- shift(tmp, shift=0:(width(peak)-1)) mcols(tmp) <- peak12725 colnames(mcols(tmp)) <- "score" peak12725 <- new("track", dat=tmp, name="peak12725", type="data", format="BED") } trs <- geneModelFromTxdb(TxDb.Hsapiens.UCSC.hg19.knownGene, org.Hs.eg.db, gr) names(trs) <- paste(sapply(trs, function(.ele) .ele@name), names(trs), sep=":") optSty <- optimizeStyle(trackList(peak12725, trs, heightDist = c(.3, .7)), theme="bw") viewTracks(optSty$tracks, gr=gr, viewerStyle=optSty$style)