annotate region_motif_intersect.r @ 36:b868314320e2 draft

Uploaded
author jeremyjliu
date Fri, 06 Mar 2015 13:24:03 -0500
parents 4ce22698acb0
children 76e108a19e64
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
21
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
1 # Name: region_motif_intersect.r
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
2 # Description: Takes a bed file of target regions and counts intersections
31
0e031bf8c5b1 Uploaded
jeremyjliu
parents: 21
diff changeset
3 # of each motif (in separately installed tabix database) and target regions.
35
4ce22698acb0 Uploaded
jeremyjliu
parents: 31
diff changeset
4 # Author: Jeremy Liu
21
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
5 # Email: jeremy.liu@yale.edu
35
4ce22698acb0 Uploaded
jeremyjliu
parents: 31
diff changeset
6 # Date: 15/02/11
21
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
7 # Note: This script is meant to be invoked with the following command
31
0e031bf8c5b1 Uploaded
jeremyjliu
parents: 21
diff changeset
8 # R --slave --vanilla -f ./region_motif_intersect.r --args <db_bgz> <db_tbi> <inbed> <outtab>
0e031bf8c5b1 Uploaded
jeremyjliu
parents: 21
diff changeset
9 # Dependencies: region_motif_data_manager
21
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
10
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
11 # Auxiliary function to concatenate multiple strings
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
12 concat <- function(...) {
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
13 input_list <- list(...)
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
14 return(paste(input_list, sep="", collapse=""))
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
15 }
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
16
31
0e031bf8c5b1 Uploaded
jeremyjliu
parents: 21
diff changeset
17 # Retrive motif database path
21
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
18 args <- commandArgs()
36
b868314320e2 Uploaded
jeremyjliu
parents: 35
diff changeset
19 motifDB_bgz = unlist(strsplit(args[7], ','))[1] # Handles duplicate entries in data table
b868314320e2 Uploaded
jeremyjliu
parents: 35
diff changeset
20 motifDB_tbi = unlist(strsplit(args[8], ','))[1] # Just takes the first one
21
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
21
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
22 # Set input and reference files, comment to toggle commmand line arguments
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
23 inBed = args[9]
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
24 outTab = args[10]
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
25
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
26 # Auxiliary function to read in BED file
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
27 read_bed <- function(file) {
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
28 return(read.table(file, sep="\t", stringsAsFactors=FALSE))
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
29 }
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
30
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
31 startTime = Sys.time()
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
32 cat("Running ... Started at:", format(startTime, "%a %b %d %X %Y"), "...\n")
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
33
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
34 # Load dependencies
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
35 cat("Loading dependencies...\n")
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
36 suppressPackageStartupMessages(library(Rsamtools, quietly=TRUE))
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
37
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
38 # Initializing hash table (as env) with motif names and loading tabix file
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
39 cat("Loading motif database and initializing hash table...\n")
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
40 motifTable = new.env()
31
0e031bf8c5b1 Uploaded
jeremyjliu
parents: 21
diff changeset
41 motifTbx <- TabixFile(motifDB_bgz)
21
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
42
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
43 # Loading input bed file, convert integer columns to numeric, name columns
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
44 cat("Loading region file...\n")
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
45 regionsDF = read_bed(inBed)
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
46 dfTemp = sapply(regionsDF, is.integer)
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
47 regionsDF[dfTemp] = lapply(regionsDF[dfTemp], as.numeric)
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
48 names(regionsDF)[names(regionsDF) == "V1"] = "chr"
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
49 names(regionsDF)[names(regionsDF) == "V2"] = "start"
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
50 names(regionsDF)[names(regionsDF) == "V3"] = "end"
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
51
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
52 # Filtering regions to exclude chromosomes not in motif database
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
53 cat("Determining intersection counts...\n")
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
54 motifTbxChrs = seqnamesTabix(motifTbx)
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
55 regionsDFFilter = subset(regionsDF, chr %in% motifTbxChrs)
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
56
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
57 # Loading regions into GRanges object and scanning motif tabix database
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
58 # Region end is incremented by 1 since scanTabix querying is inclusive for
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
59 # position start but exclusive for position end.
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
60 param = GRanges(regionsDFFilter$chr, IRanges(regionsDFFilter$start,
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
61 end=regionsDFFilter$end + 1))
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
62 regionsIntersects = scanTabix(motifTbx, param=param)
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
63
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
64 # Parsing result list and updating motif count hash table
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
65 cat("Parsing result list...\n")
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
66 for(regionIntersects in regionsIntersects) {
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
67 for(regionIntersect in strsplit(regionIntersects, " ")) {
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
68 intersectMotif = strsplit(regionIntersect, "\t")[[1]][4]
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
69 if(is.null(motifTable[[intersectMotif]])) {
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
70 motifTable[[intersectMotif]] = 1
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
71 } else {
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
72 motifTable[[intersectMotif]] = motifTable[[intersectMotif]] + 1
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
73 }
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
74 }
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
75 }
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
76
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
77 # Converting motif count hash table to an integer vector for output
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
78 counts = integer(length = length(ls(motifTable)))
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
79 names(counts) = ls(motifTable)
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
80 for(motifName in ls(motifTable)) {
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
81 counts[motifName] = as.integer(motifTable[[motifName]])
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
82 }
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
83
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
84 # Outputting intersection counts to tab delineated file
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
85 cat("Outputting to file...\n")
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
86 write.table(counts, outTab, quote=FALSE, sep="\t", row.names=TRUE, col.names=FALSE)
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
87 cat("Done. Job started at:", format(startTime, "%a %b %d %X %Y."),
2c909bfdd090 Uploaded
jeremyjliu
parents:
diff changeset
88 "Job ended at:", format(Sys.time(), "%a %b %d %X %Y."), "\n")