In CharlesJB/peakimport: Import narrowPeak and broadPeak files

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Vignette setup

For this vignette, we will import narrowPeak and broadPeak files. We will then show how to annotate the results with the r Biocpkg("ChIPseeker") package. To complete this demo, you will need the following packages:

library(peakimport)
require(ChIPseeker)
require(TxDb.Hsapiens.UCSC.hg19.knownGene)
require(org.Hs.eg.db)
require(purrr)

Importing files

In this section we will show how to import both narrowPeak and broadPeak files. We will first show how to import each type separately, and then we will show how we can use lapply or map to import multiple files with a single command.

narrowPeak files

The narrowPeak file format is described on the UCSC webpage.

It is not a valid format and cannot be directly imported using the import function from rtracklayer.

narrowPeak_file <- get_demo_file("narrowPeak")
gr_narrowPeak <- import_peaks(narrowPeak_file)
gr_narrowPeak

The results is a GRanges object.

broadPeak files

The broadPeak file format is described on the UCSC webpage.

Similar to the narrowPeak format, it is not a valid bed file and will need to be imported accordingly.

broadPeak_file <- get_demo_file("broadPeak")
gr_broadPeak <- import_peaks(broadPeak_file)
gr_broadPeak

The results is also a GRanges object.

Importing multiple files

We can use functions such as lapply or map (from the r CRANpkg("purrr") package) to import multiple files in a single command:

multiple_files <- c(narrowPeak_file, broadPeak_file)
names(multiple_files) <- c("narrowPeak", "broadPeak")
grl <- lapply(multiple_files, import_peaks)
# or
grl <- map(multiple_files, import_peaks)
grl

Both lapply and map will return the same result: a list of GRanges.

The import_peaks function will determine if the file is a bed, narrowPeak or broadPeak file based on the file extension and will select the correct import strategy accordingly. It is important for the file extensions to be exactly .narrowPeak and .broadPeak.

Annotating results

We will use the r Biocpkg("ChIPseeker") package to annotate the results from import_peaks:

txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

anno_narrowPeak <- annotatePeak(gr_narrowPeak,
                TxDb = txdb,
                annoDb = "org.Hs.eg.db")
anno_narrowPeak
anno_broadPeak <- annotatePeak(gr_broadPeak,
                TxDb = txdb,
                annoDb = "org.Hs.eg.db")
anno_broadPeak

The result can be converted in the GRanges format with the as.GRanges function:

gr_anno_narrowPeak <- as.GRanges(anno_narrowPeak)
gr_anno_narrowPeak
gr_anno_broadPeak <- as.GRanges(anno_broadPeak)
gr_anno_broadPeak

Finally, we can annotate multiple peaks in a single command with lapply or map:

anno_grl <- map(grl, annotatePeak, TxDb = txdb, annoDb = "org.Hs.eg.db")
gr_anno_grl <- map(anno_grl, as.GRanges)
gr_anno_grl

sessionInfo()

sessionInfo()

CharlesJB/peakimport documentation built on June 3, 2019, 5:47 p.m.