annotate batch-consistency-plot.r @ 19:11269f3b68a0 draft

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