inst/doc/quickStart.R

In ChIPpeakAnno: Batch annotation of the peaks identified from either ChIP-seq, ChIP-chip experiments or any experiments resulted in large number of chromosome ranges

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