view annotate_peak/annotatePeak.R @ 15:53df7871db21 draft

Uploaded
author testtool
date Thu, 16 Mar 2017 11:03:24 -0400
parents 4a489b0d247f
children ee3fc24b2108
line wrap: on
line source

options(warn = -1)
args <- commandArgs(trailingOnly = TRUE)
DMR = args[1]
annoPeakTable = args[2]

#DMR <- ("test-data/DMR.bed")
DMRInfo = read.table(
  DMR,
  header = FALSE,
  sep = "\t",
  stringsAsFactors = FALSE,
  quote = ""
)

peaks <- GRanges(seqnames = DMRInfo[, 1],
                    ranges = IRanges
                    (start = DMRInfo[, 2], end = DMRInfo[, 3]))
peaks[1:2]
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("EnsDb.Hsapiens.v75")
library(EnsDb.Hsapiens.v75)
annoData <- toGRanges(EnsDb.Hsapiens.v75)
annoData[1:2]

## keep the seqnames in the same style
seqlevelsStyle(peaks) <- seqlevelsStyle(annoData)
## do annotation by nearest TSS
anno <- annotatePeakInBatch(peaks, AnnotationData=annoData)
anno[1:2]
#annoPeakTable <- ('test-data/ChIPPeak.csv')
write.csv(anno, annoPeakTable, row.names = FALSE)