Mercurial > repos > testtool > annotate_peak

comparison annotate_peak/annotatePeak.R @ 14:4a489b0d247f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author testtool
date Thu, 16 Mar 2017 10:59:08 -0400
parents 77c0aee75ac6
children 53df7871db21
comparison
equal deleted inserted replaced
13:77c0aee75ac6 14:4a489b0d247f
1 require("ChIPseeker", quietly = TRUE)
2 require("ChIPpeakAnno", quietly = TRUE)
3
4 options(warn = -1) 1 options(warn = -1)
5 options("download.file.method"="wget")
6
7 args <- commandArgs(trailingOnly = TRUE) 2 args <- commandArgs(trailingOnly = TRUE)
8 DMR = args[1] 3 DMR = args[1]
9 annoPeakTable = args[2] 4 annoPeakTable = args[2]
10 5
11 DMR <- ("test-data/DMR.bed") 6 DMR <- ("test-data/DMR.bed")
15 sep = "\t", 10 sep = "\t",
16 stringsAsFactors = FALSE, 11 stringsAsFactors = FALSE,
17 quote = "" 12 quote = ""
18 ) 13 )
19 14
20 DMRPeaks <- GRanges(seqnames = DMRInfo[, 1], 15 peaks <- GRanges(seqnames = DMRInfo[, 1],
21 ranges = IRanges 16 ranges = IRanges
22 (start = DMRInfo[, 2], end = DMRInfo[, 3])) 17 (start = DMRInfo[, 2], end = DMRInfo[, 3]))
18 peaks[1:2]
19 ## try http:// if https:// URLs are not supported
20 source("https://bioconductor.org/biocLite.R")
21 biocLite("EnsDb.Hsapiens.v75")
22 library(EnsDb.Hsapiens.v75)
23 annoData <- toGRanges(EnsDb.Hsapiens.v75)
24 annoData[1:2]
23 25
24 annotatePeak <- 26 ## keep the seqnames in the same style
25 as.data.frame(annotatePeak(DMRPeaks, level = "gene", annoDb = "org.Hs.eg.db")) 27 seqlevelsStyle(peaks) <- seqlevelsStyle(annoData)
28 ## do annotation by nearest TSS
29 anno <- annotatePeakInBatch(peaks, AnnotationData=annoData)
30 anno[1:2]
31 annoPeakTable <- ('test-data/ChIPPeak.csv')
32 write.csv(anno, annoPeakTable, row.names = FALSE)
26 33
27 annoPeakTable <- ('test-data/annoPeakTable.csv') 34
28 write.csv(annotatePeak, annoPeakTable, row.names = FALSE)