|
5
|
1 ##############################################################################
|
|
|
2
|
|
|
3 # Modified 06/29/12: Kar Ming Chu
|
|
|
4 # Modified to work with Galaxy
|
|
|
5
|
|
|
6 # Usage: Rscript batch-consistency-analysis.r peakfile1 peakfile2 half.width overlap.ratio is.broadpeak sig.value gtable r.output overlap.output npeaks.output em.sav.output uri.sav.output
|
|
|
7
|
|
|
8 # Changes:
|
|
|
9 # - Appended parameter for input gnome table called gtable
|
|
|
10 # - Appended parameter for specifying Rout output file name (required by Galaxy)
|
|
|
11 # - Appended parameter for specifying Peak overlap output file name (required by Galaxy)
|
|
|
12 # - Appended parameter for specifying Npeak above IDR output file name (required by Galaxy)
|
|
|
13 # - Removed parameter outfile.prefix since main output files are replaced with strict naming
|
|
|
14 # - Appended parameter for specifying em.sav output file (for use with batch-consistency-plot.r)
|
|
|
15 # - Appended parameter for specifying uri.sav output file (for use with batch-consistency-plot.r)
|
|
|
16
|
|
|
17 ##############################################################################
|
|
|
18
|
|
|
19 # modified 3-29-10: Qunhua Li
|
|
|
20 # add 2 columns in the output of "-overlapped-peaks.txt": local.idr and IDR
|
|
|
21
|
|
|
22 # 01-20-2010 Qunhua Li
|
|
|
23 #
|
|
|
24 # This program performs consistency analysis for a pair of peak calling outputs
|
|
|
25 # It takes narrowPeak or broadPeak formats.
|
|
|
26 #
|
|
|
27 # usage: Rscript batch-consistency-analysis2.r peakfile1 peakfile2 half.width outfile.prefix overlap.ratio is.broadpeak sig.value
|
|
|
28 #
|
|
|
29 # peakfile1 and peakfile2 : the output from peak callers in narrowPeak or broadPeak format
|
|
|
30 # half.width: -1 if using the reported peak width,
|
|
|
31 # a numerical value to truncate the peaks to
|
|
|
32 # outfile.prefix: prefix of output file
|
|
|
33 # overlap.ratio: a value between 0 and 1. It controls how much overlaps two peaks need to have to be called as calling the same region. It is the ratio of overlap / short peak of the two. When setting at 0, it means as long as overlapped width >=1bp, two peaks are deemed as calling the same region.
|
|
|
34 # is.broadpeak: a logical value. If broadpeak is used, set as T; if narrowpeak is used, set as F
|
|
|
35 # sig.value: type of significant values, "q.value", "p.value" or "signal.value" (default, i.e. fold of enrichment)
|
|
|
36
|
|
|
37 args <- commandArgs(trailingOnly=T)
|
|
|
38
|
|
|
39 # consistency between peakfile1 and peakfile2
|
|
|
40 #input1.dir <- args[1]
|
|
|
41 #input2.dir <- args[2] # directories of the two input files
|
|
|
42 peakfile1 <- args[1]
|
|
|
43 peakfile2 <- args[2]
|
|
|
44
|
|
|
45 if(as.numeric(args[3])==-1){ # enter -1 when using the reported length
|
|
|
46 half.width <- NULL
|
|
|
47 }else{
|
|
|
48 half.width <- as.numeric(args[3])
|
|
|
49 }
|
|
|
50
|
|
|
51 overlap.ratio <- args[4]
|
|
|
52
|
|
|
53 if(args[5] == "T"){
|
|
|
54 is.broadpeak <- T
|
|
|
55 }else{
|
|
|
56 is.broadpeak <- F
|
|
|
57 }
|
|
|
58
|
|
|
59 sig.value <- args[6]
|
|
|
60
|
|
|
61
|
|
|
62 #dir1 <- "~/ENCODE/anshul/data/"
|
|
|
63 #dir2 <- dir1
|
|
|
64 #peakfile1 <- "../data/SPP.YaleRep1Gm12878Cfos.VS.Gm12878Input.PointPeak.narrowPeak"
|
|
|
65 #peakfile2 <- "../data/SPP.YaleRep3Gm12878Cfos.VS.Gm12878Input.PointPeak.narrowPeak"
|
|
|
66 #half.width <- NULL
|
|
|
67 #overlap.ratio <- 0.1
|
|
|
68 #sig.value <- "signal.value"
|
|
|
69
|
|
|
70
|
|
|
71 source("/mnt/galaxyTools/galaxy-central/tools/modENCODE_DCC_tools/idr/functions-all-clayton-12-13.r")
|
|
|
72
|
|
|
73 # read the length of the chromosomes, which will be used to concatenate chr's
|
|
|
74 # chr.file <- "genome_table.txt"
|
|
|
75 # args[7] is the gtable
|
|
|
76 chr.file <- args[7]
|
|
|
77
|
|
|
78 chr.size <- read.table(chr.file)
|
|
|
79
|
|
|
80 # setting output files
|
|
|
81 r.output <- args[8]
|
|
|
82 overlap.output <- args[9]
|
|
|
83 npeaks.output <- args[10]
|
|
|
84 em.sav.output <- args[11]
|
|
|
85 uri.sav.output <- args[12]
|
|
|
86
|
|
|
87 # sink(paste(output.prefix, "-Rout.txt", sep=""))
|
|
|
88 sink(r.output)
|
|
|
89
|
|
|
90 ############# process the data
|
|
|
91 cat("is.broadpeak", is.broadpeak, "\n")
|
|
|
92 # process data, summit: the representation of the location of summit
|
|
|
93 rep1 <- process.narrowpeak(paste(peakfile1, sep=""), chr.size, half.width=half.width, summit="offset", broadpeak=is.broadpeak)
|
|
|
94 rep2 <- process.narrowpeak(paste(peakfile2, sep=""), chr.size, half.width=half.width, summit="offset", broadpeak=is.broadpeak)
|
|
|
95
|
|
|
96 cat(paste("read", peakfile1, ": ", nrow(rep1$data.ori), "peaks\n", nrow(rep1$data.cleaned), "peaks are left after cleaning\n", peakfile2, ": ", nrow(rep2$data.ori), "peaks\n", nrow(rep2$data.cleaned), " peaks are left after cleaning"))
|
|
|
97
|
|
|
98 if(args[3]==-1){
|
|
|
99 cat(paste("half.width=", "reported", "\n"))
|
|
|
100 }else{
|
|
|
101 cat(paste("half.width=", half.width, "\n"))
|
|
|
102 }
|
|
|
103 cat(paste("significant measure=", sig.value, "\n"))
|
|
|
104
|
|
|
105 # compute correspondence profile (URI)
|
|
|
106 uri.output <- compute.pair.uri(rep1$data.cleaned, rep2$data.cleaned, sig.value1=sig.value, sig.value2=sig.value, overlap.ratio=overlap.ratio)
|
|
|
107
|
|
|
108 #uri.output <- compute.pair.uri(rep1$data.cleaned, rep2$data.cleaned)
|
|
|
109
|
|
|
110 cat(paste("URI is done\n"))
|
|
|
111
|
|
|
112 # save output
|
|
|
113 # save(uri.output, file=paste(output.prefix, "-uri.sav", sep=""))
|
|
|
114 save(uri.output, file=uri.sav.output)
|
|
|
115 cat(paste("URI is saved at: ", uri.sav.output))
|
|
|
116
|
|
|
117
|
|
|
118 # EM procedure for inference
|
|
|
119 em.output <- fit.em(uri.output$data12.enrich, fix.rho2=T)
|
|
|
120
|
|
|
121 #em.output <- fit.2copula.em(uri.output$data12.enrich, fix.rho2=T, "gaussian")
|
|
|
122
|
|
|
123 cat(paste("EM is done\n\n"))
|
|
|
124
|
|
|
125 save(em.output, file=em.sav.output)
|
|
|
126 cat(paste("EM is saved at: ", em.sav.output))
|
|
|
127
|
|
|
128
|
|
|
129 # write em output into a file
|
|
|
130
|
|
|
131 cat(paste("EM estimation for the following files\n", peakfile1, "\n", peakfile2, "\n", sep=""))
|
|
|
132
|
|
|
133 print(em.output$em.fit$para)
|
|
|
134
|
|
|
135 # add on 3-29-10
|
|
|
136 # output both local idr and IDR
|
|
|
137 idr.local <- 1-em.output$em.fit$e.z
|
|
|
138 IDR <- c()
|
|
|
139 o <- order(idr.local)
|
|
|
140 IDR[o] <- cumsum(idr.local[o])/c(1:length(o))
|
|
|
141
|
|
|
142
|
|
|
143 write.out.data <- data.frame(chr1=em.output$data.pruned$sample1[, "chr"],
|
|
|
144 start1=em.output$data.pruned$sample1[, "start.ori"],
|
|
|
145 stop1=em.output$data.pruned$sample1[, "stop.ori"],
|
|
|
146 sig.value1=em.output$data.pruned$sample1[, "sig.value"],
|
|
|
147 chr2=em.output$data.pruned$sample2[, "chr"],
|
|
|
148 start2=em.output$data.pruned$sample2[, "start.ori"],
|
|
|
149 stop2=em.output$data.pruned$sample2[, "stop.ori"],
|
|
|
150 sig.value2=em.output$data.pruned$sample2[, "sig.value"],
|
|
|
151 idr.local=1-em.output$em.fit$e.z, IDR=IDR)
|
|
|
152
|
|
|
153 # write.table(write.out.data, file=paste(output.prefix, "-overlapped-peaks.txt", sep=""))
|
|
|
154 write.table(write.out.data, file=overlap.output)
|
|
|
155 cat(paste("Write overlapped peaks and local idr to: ", overlap.output, sep=""))
|
|
|
156
|
|
|
157 # number of peaks passing IDR range (0.01-0.25)
|
|
|
158 IDR.cutoff <- seq(0.01, 0.25, by=0.01)
|
|
|
159 idr.o <- order(write.out.data$idr.local)
|
|
|
160 idr.ordered <- write.out.data$idr.local[idr.o]
|
|
|
161 IDR.sum <- cumsum(idr.ordered)/c(1:length(idr.ordered))
|
|
|
162
|
|
|
163 IDR.count <- c()
|
|
|
164 n.cutoff <- length(IDR.cutoff)
|
|
|
165 for(i in 1:n.cutoff){
|
|
|
166 IDR.count[i] <- sum(IDR.sum <= IDR.cutoff[i])
|
|
|
167 }
|
|
|
168
|
|
|
169
|
|
|
170 # write the number of peaks passing various IDR range into a file
|
|
|
171 idr.cut <- data.frame(peakfile1, peakfile2, IDR.cutoff=IDR.cutoff, IDR.count=IDR.count)
|
|
|
172 write.table(idr.cut, file=npeaks.output, append=T, quote=F, row.names=F, col.names=F)
|
|
|
173 cat(paste("Write number of peaks above IDR cutoff [0.01, 0.25]: ","npeaks-aboveIDR.txt\n", sep=""))
|
|
|
174
|
|
|
175 mar.mean <- get.mar.mean(em.output$em.fit)
|
|
|
176
|
|
|
177 cat(paste("Marginal mean of two components:\n"))
|
|
|
178 print(mar.mean)
|
|
|
179
|
|
|
180 sink()
|
|
|
181
|
|
|
182
|
