Mercurial > repos > rnateam > chipseeker

comparison chipseeker.R @ 3:535321abf9a4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/chipseeker commit 861db0d1f76bb320f49c2501f4e656cf88d389ce
author rnateam
date Fri, 01 Jun 2018 02:13:33 -0400
parents 95f779f4adb7
children 90fe78a19b55
comparison
equal deleted inserted replaced
2:cb133602cd9b 3:535321abf9a4
4 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") 4 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
5 5
6 suppressPackageStartupMessages({ 6 suppressPackageStartupMessages({
7 library(ChIPseeker) 7 library(ChIPseeker)
8 library(GenomicFeatures) 8 library(GenomicFeatures)
9 library(rtracklayer)
9 library(optparse) 10 library(optparse)
10 }) 11 })
11 12
12 option_list <- list( 13 option_list <- list(
13 make_option(c("-i","--infile"), type="character", help="Peaks file to be annotated"), 14 make_option(c("-i","--infile"), type="character", help="Peaks file to be annotated"),
15 make_option(c("-u","--upstream"), type="integer", help="TSS upstream region"), 16 make_option(c("-u","--upstream"), type="integer", help="TSS upstream region"),
16 make_option(c("-d","--downstream"), type="integer", help="TSS downstream region"), 17 make_option(c("-d","--downstream"), type="integer", help="TSS downstream region"),
17 make_option(c("-F","--flankgeneinfo"), type="logical", help="Add flanking gene info"), 18 make_option(c("-F","--flankgeneinfo"), type="logical", help="Add flanking gene info"),
18 make_option(c("-D","--flankgenedist"), type="integer", help="Flanking gene distance"), 19 make_option(c("-D","--flankgenedist"), type="integer", help="Flanking gene distance"),
19 make_option(c("-f","--format"), type="character", help="Output format (interval or tabular)."), 20 make_option(c("-f","--format"), type="character", help="Output format (interval or tabular)."),
20 make_option(c("-p","--plots"), type="character", help="PDF of plots.") 21 make_option(c("-p","--plots"), type="logical", help="PDF of plots."),
22 make_option(c("-r","--rdata"), type="logical", help="Output RData file.")
21 ) 23 )
22 24
23 parser <- OptionParser(usage = "%prog [options] file", option_list=option_list) 25 parser <- OptionParser(usage = "%prog [options] file", option_list=option_list)
24 args = parse_args(parser) 26 args = parse_args(parser)
25 27
26 peaks = args$infile 28 peaks = args$infile
27 gtf = args$gtf 29 gtf = args$gtf
28 up = args$upstream 30 up = args$upstream
29 down = args$downstream 31 down = args$downstream
30 format = args$format 32 format = args$format
31 plots = args$plots
32 33
33 peaks <- readPeakFile(peaks) 34 peaks <- readPeakFile(peaks)
34 35
35 # Make TxDb from GTF 36 # Make TxDb from GTF
36 txdb <- makeTxDbFromGFF(gtf, format="gtf") 37 txdb <- makeTxDbFromGFF(gtf, format="gtf")
38 peakAnno <- annotatePeak(peaks, TxDb=txdb, tssRegion=c(-up, down), addFlankGeneInfo=args$flankgeneinfo, flankDistance=args$flankgenedist) 39 peakAnno <- annotatePeak(peaks, TxDb=txdb, tssRegion=c(-up, down), addFlankGeneInfo=args$flankgeneinfo, flankDistance=args$flankgenedist)
39 } else { 40 } else {
40 peakAnno <- annotatePeak(peaks, TxDb=txdb, tssRegion=c(-up, down)) 41 peakAnno <- annotatePeak(peaks, TxDb=txdb, tssRegion=c(-up, down))
41 } 42 }
42 43
44 # Add gene name
45 features <- import(gtf, format="gtf")
46 ann <- unique(mcols(features)[, c("gene_id", "gene_name")])
47 res <- as.data.frame(peakAnno)
48 res <- merge(res, ann, by.x="geneId", by.y="gene_id")
49 names(res)[names(res) == "gene_name"] <- "geneName"
50
51 #Extract metadata cols, 1st is geneId, rest should be from col 7 to end
52 metacols <- res[, c(7:ncol(res), 1)]
43 # Convert from 1-based to 0-based format 53 # Convert from 1-based to 0-based format
44 res <- as.GRanges(peakAnno)
45 metacols <- mcols(res)
46 if (format == "interval") { 54 if (format == "interval") {
47 metacols <- apply(as.data.frame(metacols), 1, function(col) paste(col, collapse="|")) 55 metacols <- apply(as.data.frame(metacols), 1, function(col) paste(col, collapse="|"))
48 resout <- data.frame(Chrom=seqnames(res), 56 resout <- data.frame(Chrom=res$seqnames,
49 Start=start(res) - 1, 57 Start=res$start - 1,
50 End=end(res), 58 End=res$end,
51 Comment=metacols) 59 Comment=metacols)
52 } else { 60 } else {
53 resout <- data.frame(Chrom=seqnames(res), 61 resout <- data.frame(Chrom=res$seqnames,
54 Start=start(res) - 1, 62 Start=res$start - 1,
55 End=end(res), 63 End=res$end,
56 metacols) 64 metacols)
57 } 65 }
58 66
59 write.table(resout, file="out.tab", sep="\t", row.names=FALSE, quote=FALSE) 67 write.table(resout, file="out.tab", sep="\t", row.names=FALSE, quote=FALSE)
60 68
61 if (!is.null(plots)) { 69 if (!is.null(args$plots)) {
62 pdf("out.pdf", width=14) 70 pdf("out.pdf", width=14)
63 plotAnnoPie(peakAnno) 71 plotAnnoPie(peakAnno)
64 plotAnnoBar(peakAnno) 72 plotAnnoBar(peakAnno)
65 vennpie(peakAnno) 73 vennpie(peakAnno)
66 upsetplot(peakAnno) 74 upsetplot(peakAnno)
67 plotDistToTSS(peakAnno, title="Distribution of transcription factor-binding loci\nrelative to TSS") 75 plotDistToTSS(peakAnno, title="Distribution of transcription factor-binding loci\nrelative to TSS")
68 dev.off() 76 dev.off()
69 } 77 }
78
79 ## Output RData file
80
81 if (!is.null(args$rdata)) {
82 save.image(file = "ChIPseeker_analysis.RData")
83 }