Extract ChIP Seq profiles using metagene

With metagene package, it is easy to extract the genomic profiles from BAM files.

It can be done in few simple steps:

1. Loading metagene package

The metagene package is mandatory. So, metagene should be installed and loaded.

library(metagene)

2. Preparing BAM files

Aligned sequences are usually stored in BAM files. The BAM file(s) corresponding of the ChIP-Seq of interest have to be loaded.

bam_files <- c(system.file("extdata/align1_rep1.bam", package = "metagene"),
               system.file("extdata/align1_rep2.bam", package = "metagene"))

2. Selecting regions of interest

Usually, ChIP-Seq analysis are only done in regions of interest, such as TSS. So, regions of interest can be specified so that only those regions are extracted from BAM files. In this example, regions are kept in a text file and tranformed into GRanges objects. However, regions can come under different formats, such as BED files and GRanges objects.

library(GenomicRanges)
regions <- system.file("extdata/list1.BED", package = "metagene")
regions <- read.table(regions, header = FALSE)
names(regions)<-c("chromosome", "chromStart", "chromEnd")
regions <- makeGRangesFromDataFrame(regions, start.field = "chromStart",
               end.field = "chromEnd")

3. Creating metagene object

A metagene object is created using the BAM files and regions as parameters.

mg <- metagene$new(regions = regions, bam_files = bam_files, cores = 4)

4. Extracting ChIP-Seq profiles

The metagene object contains variable which allow acces to coverage information.

profiles <- mg$coverages

The profiles is a SimpleRleList object.

5. Converting into ChIP-Seq profiles vector

The profiles object can easily be converted into coverage vectors.

bamfile_align1_rep1 <- bam_files[1]
profiles_align1_rep1 <- profiles[[bamfile_align1_rep1]][[1]]
profiles_align1_rep1 <- lapply(profiles_align1_rep1, as.numeric)

bamfile_align1_rep2 <- bam_files[2]
profiles_align1_rep2 <- profiles[[bamfile_align1_rep2]][[1]]
profiles_align1_rep2 <- lapply(profiles_align1_rep2, as.numeric)