Mercurial > repos > jeremyjliu > region_motif_enrichment

comparison region_motif_compare.r @ 3:cab2db9d058b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author jeremyjliu
date Sat, 16 May 2015 22:35:26 -0400
parents
children 4803f5186f1a
comparison
equal deleted inserted replaced
2:2538677b6004 3:cab2db9d058b
1 # Name: region_motif_compare.r
2 # Description: Reads in two count files and determines enriched and depleted
3 # motifs (or any location based feature) based on poisson tests and gc
4 # corrections. All enrichment ratios relative to overall count / gc ratios.
5 # Author: Jeremy liu
6 # Email: jeremy.liu@yale.edu
7 # Date: 15/02/11
8 # Note: This script can be invoked with the following command
9 # R --slave --vanilla -f ./region_motif_compare.r --args <workingdir> <pwm_file>
10 # <intab1> <intab2> <enriched_tab> <depleted_tab> <plots_png>
11 # <workingdir> is the directory where plotting.r is saved
12 # Dependencies: region_motif_data_manager, plotting.r,
13
14 # Auxiliary function to concatenate multiple strings
15 concat <- function(...) {
16 input_list <- list(...)
17 return(paste(input_list, sep="", collapse=""))
18 }
19
20 # Supress all warning messages to prevent Galaxy treating warnings as errors
21 options(warn=-1)
22
23 # Set common and data directories
24 args <- commandArgs()
25 workingDir = args[7]
26 pwmFile = unlist(strsplit(args[8], ','))[1] # If duplicate entires, take first one
27
28 # Set input and reference files
29 inTab1 = args[9]
30 inTab2 = args[10]
31 enrichTab = args[11]
32 depleteTab = args[12]
33 plotsPng = args[13]
34
35 # Load dependencies
36 source(concat(workingDir, "/plotting.r"))
37
38 # Auxiliary function to read in tab file and prepare the data
39 read_tsv <- function(file) {
40 data = read.table(file, sep="\t", stringsAsFactors=FALSE)
41 names(data)[names(data) == "V1"] = "motif"
42 names(data)[names(data) == "V2"] = "counts"
43 return(data)
44 }
45
46 startTime = Sys.time()
47 cat("Running ... Started at:", format(startTime, "%a %b %d %X %Y"), "...\n")
48
49 # Loading motif position weight matrix (pwm) file
50 cat("Loading motif postion weight matrices...\n")
51 lines = scan(pwmFile, what="character", sep="\n", quiet=TRUE)
52 indices = which(grepl("MOTIF", lines))
53 names(indices) = lapply(indices, function(i) {
54 nameline = lines[i]
55 name = substr(nameline, 7, nchar(nameline))
56 })
57
58 pwms = sapply(indices, function(i) {
59 infoline = unlist(strsplit(lines[i+1], " "))
60 alength = as.numeric(infoline[4])
61 width = as.numeric(infoline[6])
62 subset = lines[(i+2):(i+2+width-1)]
63 motiflines = strsplit(subset, " ")
64 motif = t(do.call(rbind, motiflines))
65 motif = apply(motif, 2, as.numeric)
66 }, simplify=FALSE, USE.NAMES=TRUE)
67
68 # Loading input tab files
69 cat("Loading and reading input region motif count files...\n")
70 region1DF = read_tsv(inTab1)
71 region2DF = read_tsv(inTab2)
72 region1Counts = region1DF$counts
73 region2Counts = region2DF$counts
74 names(region1Counts) = region1DF$motif
75 names(region2Counts) = region2DF$motif
76
77 # Processing count vectors to account for missing 0 count motifs, then sorting
78 cat("Performing 0 count correction and sorting...\n")
79 allNames = union(names(region1Counts), names(region2Counts))
80 region1Diff = setdiff(allNames, names(region1Counts))
81 region2Diff = setdiff(allNames, names(region2Counts))
82 addCounts1 = rep(0, length(region1Diff))
83 addCounts2 = rep(0, length(region2Diff))
84 names(addCounts1) = region1Diff
85 names(addCounts2) = region2Diff
86 newCounts1 = append(region1Counts, addCounts1)
87 newCounts2 = append(region2Counts, addCounts2)
88 region1Counts = newCounts1[sort.int(names(newCounts1), index.return=TRUE)$ix]
89 region2Counts = newCounts2[sort.int(names(newCounts2), index.return=TRUE)$ix]
90
91 # Generate gc content matrix
92 gc = sapply(pwms, function(i) mean(i[2:3,3:18]))
93
94 # Apply poisson test, calculate p and q values, and filter significant results
95 cat("Applying poisson test...\n")
96 rValue = sum(region2Counts) / sum(region1Counts)
97 pValue = sapply(seq(along=region1Counts), function(i) {
98 poisson.test(c(region1Counts[i], region2Counts[i]), r=1/rValue)$p.value
99 })
100 qValue = p.adjust(pValue, "fdr")
101 indices = which(qValue<0.1 & abs(log2(region1Counts/region2Counts/rValue))>log2(1.5))
102
103 # Setting up output diagnostic plots, 4 in 1 png image
104 png(plotsPng, width=800, height=800)
105 xlab = "region1_count"
106 ylab = "region2_count"
107 lim = c(0.5, 5000)
108 layout(matrix(1:4, ncol=2))
109 par(mar=c(5, 5, 5, 1))
110
111 # Plot all motif counts along the linear correlation coefficient
112 plot.scatter(region1Counts+0.5, region2Counts+0.5, log="xy", xlab=xlab, ylab=ylab,
113 cex.lab=2.2, cex.axis=1.8, xlim=lim, ylim=lim*rValue)
114 abline(0, rValue, untf=T)
115 abline(0, rValue*2, untf=T, lty=2)
116 abline(0, rValue/2, untf=T, lty=2)
117
118 # Plot enriched and depleted motifs in red, housed in second plot
119 plot.scatter(region1Counts+0.5, region2Counts+0.5, log="xy", xlab=xlab, ylab=ylab,
120 cex.lab=2.2, cex.axis=1.8, xlim=lim, ylim=lim*rValue)
121 points(region1Counts[indices]+0.5, region2Counts[indices]+0.5, col="red")
122 abline(0, rValue, untf=T)
123 abline(0, rValue*2, untf=T, lty=2)
124 abline(0, rValue/2, untf=T, lty=2)
125
126 # Apply and plot gc correction and loess curve
127 cat("Applying gc correction, rerunning poisson test...\n")
128 ind = which(region1Counts>5)
129 gc = gc[names(region2Counts)] # Reorder the indices of pwms to match input data
130 lo = plot.scatter(gc,log2(region2Counts/region1Counts),draw.loess=T,
131 xlab="gc content of motif",ylab=paste("log2(",ylab,"/",xlab,")"),
132 cex.lab=2.2,cex.axis=1.8,ind=ind) # This function is in plotting.r
133 gcCorrection = 2^approx(lo$loess,xout=gc,rule=2)$y
134
135 # Recalculate p and q values, and filter for significant entries
136 pValueGC = sapply(seq(along=region1Counts),function(i) {
137 poisson.test(c(region1Counts[i],region2Counts[i]),r=1/gcCorrection[i])$p.value
138 })
139 qValueGC=p.adjust(pValueGC,"fdr")
140 indicesGC = which(qValueGC<0.1 & abs(log2(region1Counts/region2Counts*gcCorrection))>log2(1.5))
141
142 # Plot gc corrected motif counts
143 plot.scatter(region1Counts+0.5, (region2Counts+0.5)/gcCorrection, log="xy",
144 xlab=xlab, ylab=paste(ylab,"(normalized)"), cex.lab=2.2, cex.axis=1.8,
145 xlim=lim, ylim=lim)
146 points(region1Counts[indicesGC]+0.5,
147 (region2Counts[indicesGC]+0.5)/gcCorrection[indicesGC], col="red")
148 abline(0,1)
149 abline(0,1*2,untf=T,lty=2)
150 abline(0,1/2,untf=T,lty=2)
151
152 # Trim results, compile statistics and output to file
153 # Only does so if significant results are computed
154 if(length(indicesGC) > 0) {
155 # Calculate expected counts and enrichment ratios
156 cat("Calculating statistics...\n")
157 nullExpect = region1Counts * gcCorrection
158 enrichment = region2Counts / nullExpect
159
160 # Reorder selected indices in ascending pvalue
161 cat("Reordering by ascending pvalue...\n")
162 indicesReorder = indicesGC[order(pValueGC[indicesGC])]
163
164 # Combine data into one data frame and output to two files
165 cat("Splitting and outputting data...\n")
166 outDF = data.frame(motif=names(pValueGC), p=as.numeric(pValueGC), q=qValueGC,
167 stringsAsFactors=F, region_1_count=region1Counts,
168 null_expectation=round(nullExpect,2), region_2_count=region2Counts,
169 enrichment=enrichment)[indicesReorder,]
170 names(outDF)[which(names(outDF)=="region_1_count")]=xlab
171 names(outDF)[which(names(outDF)=="region_2_count")]=ylab
172 indicesEnrich = which(outDF$enrichment>1)
173 indicesDeplete = which(outDF$enrichment<1)
174 outDF$enrichment = ifelse(outDF$enrichment>1,
175 round(outDF$enrichment,3),
176 paste("1/",round(1/outDF$enrichment,3)))
177 write.table(outDF[indicesEnrich,], file=enrichTab, quote=FALSE,
178 sep="\t", append=FALSE, row.names=FALSE, col.names=TRUE)
179 write.table(outDF[indicesDeplete,], file=depleteTab, quote=FALSE,
180 sep="\t", append=FALSE, row.names=FALSE, col.names=TRUE)
181 }
182
183 # Catch display messages and output timing information
184 catchMessage = dev.off()
185 cat("Done. Job started at:", format(startTime, "%a %b %d %X %Y."),
186 "Job ended at:", format(Sys.time(), "%a %b %d %X %Y."), "\n")