|
4
|
1 # 1-20-10 Qunhua Li
|
|
|
2 #
|
|
|
3 # This program first plots correspondence curve and IDR threshold plot
|
|
|
4 # (i.e. number of selected peaks vs IDR) for each pair of sample
|
|
|
5 #
|
|
|
6 # usage:
|
|
|
7 # Rscript batch-consistency-plot-merged.r [npairs] [output.dir] [input.file.prefix 1, 2, 3 ...]
|
|
|
8 # [npairs]: integer, number of consistency analyses
|
|
|
9 # (e.g. if 2 replicates, npairs=1, if 3 replicates, npairs=3
|
|
|
10 # [output.prefix]: output directory and file name prefix for plot eg. /plots/idrPlot
|
|
|
11 # [input.file.prefix 1, 2, 3]: prefix for the output from batch-consistency-analysis2. They are the input files for merged analysis see below for examples (i.e. saved.file.prefix). It can be multiple files
|
|
|
12 #
|
|
|
13
|
|
|
14 args <- commandArgs(trailingOnly=T)
|
|
|
15 npair <- args[1] # number of curves to plot on the same figure
|
|
|
16 output.file.prefix <- args[2] # file name for plot, generated from script at the outer level
|
|
|
17 df.txt <- 10
|
|
|
18 ntemp <- as.numeric(npair)
|
|
|
19 saved.file.prefix <- list() # identifier of filenames that contain the em and URI results
|
|
|
20 source("/mnt/galaxyTools/galaxy-central/tools/modENCODE_DCC_tools/idr/functions-all-clayton-12-13.r")
|
|
|
21
|
|
|
22 uri.list <- list()
|
|
|
23 uri.list.match <- list()
|
|
|
24 ez.list <- list()
|
|
|
25 legend.txt <- c()
|
|
|
26 em.output.list <- list()
|
|
|
27 uri.output.list <- list()
|
|
|
28
|
|
|
29 for(i in 1:npair){
|
|
|
30 saved.file.prefix[i] <- args[2+i]
|
|
|
31
|
|
|
32 load(paste(saved.file.prefix[i], "-uri.sav", sep=""))
|
|
|
33 load(paste(saved.file.prefix[i], "-em.sav", sep=""))
|
|
|
34
|
|
|
35 uri.output.list[[i]] <- uri.output
|
|
|
36 em.output.list[[i]] <- em.output
|
|
|
37
|
|
|
38 ez.list[[i]] <- get.ez.tt.all(em.output, uri.output.list[[i]]$data12.enrich$merge1,
|
|
|
39 uri.output.list[[i]]$data12.enrich$merge2) # reverse =T for error rate
|
|
|
40
|
|
|
41 # URI for all peaks
|
|
|
42 uri.list[[i]] <- uri.output$uri.n
|
|
|
43 # URI for matched peaks
|
|
|
44 uri.match <- get.uri.matched(em.output$data.pruned, df=df.txt)
|
|
|
45 uri.list.match[[i]] <- uri.match$uri.n
|
|
|
46
|
|
|
47 file.name <- unlist(strsplit(as.character(saved.file.prefix[i]), "/"))
|
|
|
48
|
|
|
49 legend.txt[i] <- paste(i, "=", file.name[length(file.name)])
|
|
|
50
|
|
|
51 }
|
|
|
52
|
|
|
53 plot.uri.file <- paste(output.file.prefix, "-plot.ps", sep="")
|
|
|
54
|
|
|
55 ############# plot and report output
|
|
|
56 # plot correspondence curve for each pair,
|
|
|
57 # plot number of selected peaks vs IDR
|
|
|
58 # plot all into 1 file
|
|
|
59 postscript(paste(output.file.prefix, "-plot.ps", sep=""))
|
|
|
60 par(mfcol=c(2,3), mar=c(5,6,4,2)+0.1)
|
|
|
61 plot.uri.group(uri.list, NULL, file.name=NULL, c(1:npair), title.txt="all peaks")
|
|
|
62 plot.uri.group(uri.list.match, NULL, file.name=NULL, c(1:npair), title.txt="matched peaks")
|
|
|
63 plot.ez.group(ez.list, plot.dir=NULL, file.name=NULL, legend.txt=c(1:npair), y.lim=c(0, 0.6))
|
|
|
64 plot(0, 1, type="n", xlim=c(0,1), ylim=c(0,1), xlab="", ylab="", xaxt="n", yaxt="n") # legends
|
|
|
65 legend(0, 1, legend.txt, cex=0.6)
|
|
|
66
|
|
|
67 dev.off()
|
|
|
68
|
