annotate 02_read_counts_tss/countsAroundTSS.R @ 2:86c0309e6a5f draft default tip

Deleted selected files
author dktanwar
date Tue, 12 Sep 2017 14:12:35 -0400
parents a108d4a59736
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
1 ## How to execute this tool
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
2 # $Rscript my_r_tool.R --input input.csv --output output.csv
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
3
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
4 # Send R errors to stderr
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
5 options(show.error.messages = F, error = function(){cat(geterrmessage(), file = stderr()); q("no", 1, F)})
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
6
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
7 # Avoid crashing Galaxy with an UTF8 error on German LC settings
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
8 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
9
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
10 options(stringAsfactors = FALSE, useFancyQuotes = FALSE)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
11
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
12 # Check for packages requirement -----
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
13 suppressMessages(source("https://bioconductor.org/biocLite.R"))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
14
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
15 list.of.packages <- c("GenomicRanges",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
16 "Rsamtools",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
17 "rtracklayer",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
18 "parallel",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
19 "foreach",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
20 "doParallel",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
21 "GenomicFeatures",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
22 "Gviz",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
23 "BSgenome.Hsapiens.UCSC.hg19",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
24 "org.Hs.eg.db",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
25 "BSgenome.Mmusculus.UCSC.mm10",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
26 "org.Mm.eg.db",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
27 "optparse",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
28 "getopt")
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
29
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
30 new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
31 if(length(new.packages)) biocLite(new.packages)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
32
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
33 # packages required
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
34 tmp <- suppressMessages(lapply(list.of.packages, require, quietly =T,
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
35 warn.conflicts = F, character.only = TRUE))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
36 rm(tmp)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
37
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
38
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
39 # Command line options ----
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
40 # option_list <- list(
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
41 # make_option(c("-b", "--bamfile"), type="character", default=NULL,
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
42 # help="bamfile for which read counts are to be calculated around TSS", metavar="character"),
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
43 # make_option(c("-e", "--genome"), type="character", default = NULL,
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
44 # help="genome: hg19/ mm10", metavar="character"),
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
45 # make_option(c("-n", "--TSSn"), type="character", default = NULL,
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
46 # help="region upstream TSS for which the read counts are to be calculated", metavar="character"),
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
47 # make_option(c("-p", "--TSSp"), type="character", default = NULL,
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
48 # help="region downstream TSS for which the read counts are to be calculated", metavar="character"),
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
49 # make_option(c("-u", "--TxDbUCSC"), type="character", default = NULL,
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
50 # help="Transcripts database: If you already have a sqlite Transcript database for hg19/ mm10, provide it, or do not enter anything. Software will create a database to use from UCSC.", metavar="character")
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
51 # )
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
52 #
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
53 # opt_parser <- OptionParser(option_list=option_list)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
54 # opt <- parse_args(opt_parser)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
55
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
56 # Take in trailing command line arguments
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
57 args <- commandArgs(trailingOnly = TRUE)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
58 # args <- c("ENCFF027UTM.bam", "hg19", 500, 500, "transcripts_mm10_UCSC.sqlite")
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
59
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
60 # Get options using the spec as defined by the enclosed list
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
61 # Options are read from the default: commandArgs(TRUE)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
62 option_specification = matrix(c(
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
63 'input1', 'i1', 2, 'character',
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
64 'input2', 'i2', 2, 'character',
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
65 'input3', 'i3', 2, 'integer',
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
66 'input4', 'i4', 2, 'integer',
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
67 'input5', 'i5', 2, 'character',
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
68 'output', 'o', 2, 'character'
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
69 ), byrow=TRUE, ncol=4);
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
70
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
71 # Parse options
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
72 options = getopt(option_specification)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
73
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
74 if (is.null(options$input1) | is.null(options$input2)){
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
75 stop("Please provide right input parameters. See help!", call.=FALSE)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
76 }
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
77
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
78
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
79 # opt$genome <- "mm10"
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
80 # opt$TxDbUCSC <- "transcripts_mm10_UCSC.sqlite"
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
81 # opt$bamfile <- "ENCFF027UTM.bam"
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
82 # opt$TSSn <- 500
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
83 # opt$TSSp <- 500
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
84
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
85 # Transcript database ----
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
86 txdb <- NULL
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
87 if(is.null(options$input5)){
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
88 txdb <- makeTxDbFromUCSC(genome = options$input2, tablename="refGene")
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
89 saveDb(txdb, file = paste("transcripts_", options$input2, "_UCSC.sqlite", sep = ""))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
90 } else{
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
91 txdb <- loadDb(options$input5)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
92 }
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
93
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
94 txs <- transcripts(txdb, columns = c("gene_id", "tx_id"))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
95 txs <- txs[sapply(split(seq_along(txs), sapply(txs$gene_id, "[", 1)),
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
96 function(x) x[which.max(width(txs)[x])])]
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
97
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
98 geneID <- txs$"gene_id"
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
99 genes <- as.character(geneID)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
100
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
101
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
102 # Bamfiles and chromosomes
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
103 # bamfiles <- list.files(path = opt$bamdir, pattern = ".bam$", full.names = T)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
104 bamfiles <- options$input1
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
105
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
106 chrs <- scanBamHeader(bamfiles)[[1]]$targets
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
107
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
108
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
109 # TSS regions for calculationg counts ----
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
110 tss <- suppressWarnings(GRanges(seqnames(txs),
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
111 IRanges(start = ifelse(as.character(strand(txs)) == "+",
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
112 start(txs), end(txs)) - as.numeric(options$input3),
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
113 width = as.numeric(options$input3) + as.numeric(options$input4)),
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
114 strand = strand(txs), seqlengths = chrs))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
115
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
116
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
117 # Functions ----
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
118 mylapply <- function(...) {
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
119 if(require(parallel) && .Platform$OS.type == "unix")
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
120 mclapply(..., mc.cores = detectCores(), mc.preschedule = F)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
121 else
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
122 lapply(...)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
123 }
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
124
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
125 coverageForChr <- function(seqname, bamfile, chrs, fragmentSize=0, ...) {
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
126 message(", ", seqname, appendLF=FALSE)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
127 shift <- round(fragmentSize/2)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
128 param <- ScanBamParam(what=c("pos","qwidth","strand"),
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
129 which=GRanges(seqname, IRanges(1, chrs[seqname])),
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
130 flag=scanBamFlag(isUnmappedQuery=FALSE))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
131 x <- scanBam(bamfile, ..., param=param)[[1]]
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
132 coverage(IRanges(ifelse(x[["strand"]]=="+", x[["pos"]]+shift,
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
133 x[["pos"]]+x[["qwidth"]]-shift), width=1), width=chrs[seqname])
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
134 }
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
135
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
136
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
137 countsForRegions <- function(bamfiles, regions, fragmentSize=0) {
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
138 names(regions) <- seq_along(regions)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
139 chrs <- seqlengths(regions)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
140 DF <- do.call(cbind, lapply(bamfiles, function(bf) {
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
141 message("counting alignments in ",basename(bf),appendLF=FALSE)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
142 cover <- mylapply(names(chrs), coverageForChr, bamfile=bf, chrs=chrs, fragmentSize=fragmentSize)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
143 names(cover) <- names(chrs)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
144 cover <- RleList(unlist(cover), compress=FALSE)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
145 vws <- Views(cover, as(regions,"RangesList"))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
146 vs <- unlist(viewSums(vws), use.names=FALSE)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
147 vs <- vs[match(names(regions), unlist(lapply(vws,names),use.names=FALSE))]
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
148 df <- DataFrame(vs)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
149 colnames(df) <- sub(".bam","",basename(bf))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
150 message(": done")
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
151 df
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
152 }))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
153 }
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
154
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
155
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
156
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
157 # Counting reads ----
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
158 cnts <- countsForRegions(bamfiles = bamfiles, regions = tss, fragmentSize = 150)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
159 cnts <- data.frame(Genes = genes, cnts, stringsAsFactors = F, check.names = F)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
160
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
161 # Annotating genes ----
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
162 anno <- NULL
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
163 if(options$input2 == "hg19"){
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
164 anno <- select(org.Hs.eg.db, as.character(cnts$Genes), c("SYMBOL", "GENENAME"))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
165 } else if (options$input2 == "mm10"){
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
166 anno <- select(org.Mm.eg.db, as.character(cnts$Genes), c("SYMBOL", "GENENAME"))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
167 } else {
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
168 stop("Please provide correct reference genome. See help!", call.=FALSE)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
169 }
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
170
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
171 cnts <- merge(anno, cnts, by.x = "ENTREZID", by.y = "Genes")
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
172
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
173 index <- which(is.na(cnts$SYMBOL))
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
174 cnts$SYMBOL[index] <- cnts$ENTREZID[index]
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
175
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
176 write.table(cnts, file = options$output, sep = "\t", row.names = F, quote = F)
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
177
a108d4a59736 Uploaded
dktanwar
parents:
diff changeset
178 sessionInfo()