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