Annotate MACS2 peaks using ChIPpeakAnno

annotate_MACS2_peaks.R

library("ChIPpeakAnno")
library("GenomicRanges")
library("org.At.tair.db")
library("TxDb.Athaliana.BioMart.plantsmart28")
library("biomaRt")

# Annotate genomic intervals in bed format using ChIPpeakAnno 
# This script was designed for Arabidopsis, but can be easily changed for
# any other organism available through biomaRt

setwd("peaks_directory")
list.files()

annotateMacsPeaks <- function(bedFile) 
{
  # resFile <- path to resFile
  # name of the out file that contains annotated results
  
  # Read in results table
  peaks <- read.table(peakFile, header=F, sep = "\t")
  
  # Convert to GRanges
  gr <- toGRanges(bedFile, format="narrowPeak", header = FALSE)
 
  # Create GRanges object with annotations from TxDb database
  annoData <- toGRanges(TxDb.Athaliana.BioMart.plantsmart28, feature="gene")
  
  # Annotate granges with the nearest TSS
  annot <- annotatePeakInBatch(gr, 
                               AnnotationData=annoData, 
                               featureType = "TSS",
                               output="nearestLocation",
                               PeakLocForDistance = "start")
  
  # Load mart
  ensembl <- useMart(biomart = "plants_mart",
                     dataset = "athaliana_eg_gene",
                     host = "plants.ensembl.org")
  
  # Add gene information
  annot <- addGeneIDs(annot, mart = ensembl, feature_id_type = "ensembl_gene_id",
                      IDs2Add = c("entrezgene", "tair_symbol", "description"))
  
  samplename <- gsub("\\.\\w+", "", bedFile)
    
  write.table(annot, paste(samplename, "_annotated.txt", sep = ""), sep = "\t", col.names = T, 
              row.names = F, quote = F)
}

files <- list.files()[grep(".narrowPeak", list.files())]

lapply( files, function(x) annotateMacsPeaks(x) )