| Previous changeset 5:48767bec000d (2013-01-17) Next changeset 7:aeea5ac02b3b (2013-01-17) |
|
Commit message:
Uploaded |
|
added:
genome_tables/genome_table.dmel.r5.32.txt genome_tables/genome_table.human.hg18.txt genome_tables/genome_table.human.hg19.txt genome_tables/genome_table.mm9.txt genome_tables/genome_table.worm.ws220.txt |
|
removed:
batch-consistency-analysis.r batch-consistency-plot.r functions-all-clayton-12-13.r idrPlotDef.xml idrPlotWrapper.sh idrToolDef.xml |
| b |
| diff -r 48767bec000d -r 16bbe277d15d batch-consistency-analysis.r --- a/batch-consistency-analysis.r Thu Jan 17 15:45:48 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| @@ -1,182 +0,0 @@ -############################################################################## - -# Modified 06/29/12: Kar Ming Chu -# Modified to work with Galaxy - -# 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 - -# Changes: -# - Appended parameter for input gnome table called gtable -# - Appended parameter for specifying Rout output file name (required by Galaxy) -# - Appended parameter for specifying Peak overlap output file name (required by Galaxy) -# - Appended parameter for specifying Npeak above IDR output file name (required by Galaxy) -# - Removed parameter outfile.prefix since main output files are replaced with strict naming -# - Appended parameter for specifying em.sav output file (for use with batch-consistency-plot.r) -# - Appended parameter for specifying uri.sav output file (for use with batch-consistency-plot.r) - -############################################################################## - -# modified 3-29-10: Qunhua Li -# add 2 columns in the output of "-overlapped-peaks.txt": local.idr and IDR - -# 01-20-2010 Qunhua Li -# -# This program performs consistency analysis for a pair of peak calling outputs -# It takes narrowPeak or broadPeak formats. -# -# usage: Rscript batch-consistency-analysis2.r peakfile1 peakfile2 half.width outfile.prefix overlap.ratio is.broadpeak sig.value -# -# peakfile1 and peakfile2 : the output from peak callers in narrowPeak or broadPeak format -# half.width: -1 if using the reported peak width, -# a numerical value to truncate the peaks to -# outfile.prefix: prefix of output file -# 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. -# is.broadpeak: a logical value. If broadpeak is used, set as T; if narrowpeak is used, set as F -# sig.value: type of significant values, "q.value", "p.value" or "signal.value" (default, i.e. fold of enrichment) - -args <- commandArgs(trailingOnly=T) - -# consistency between peakfile1 and peakfile2 -#input1.dir <- args[1] -#input2.dir <- args[2] # directories of the two input files -peakfile1 <- args[1] -peakfile2 <- args[2] - -if(as.numeric(args[3])==-1){ # enter -1 when using the reported length - half.width <- NULL -}else{ - half.width <- as.numeric(args[3]) -} - -overlap.ratio <- args[4] - -if(args[5] == "T"){ - is.broadpeak <- T -}else{ - is.broadpeak <- F -} - -sig.value <- args[6] - - -#dir1 <- "~/ENCODE/anshul/data/" -#dir2 <- dir1 -#peakfile1 <- "../data/SPP.YaleRep1Gm12878Cfos.VS.Gm12878Input.PointPeak.narrowPeak" -#peakfile2 <- "../data/SPP.YaleRep3Gm12878Cfos.VS.Gm12878Input.PointPeak.narrowPeak" -#half.width <- NULL -#overlap.ratio <- 0.1 -#sig.value <- "signal.value" - - -source("/mnt/galaxyTools/galaxy-central/tools/modENCODE_DCC_tools/idr/functions-all-clayton-12-13.r") - -# read the length of the chromosomes, which will be used to concatenate chr's -# chr.file <- "genome_table.txt" -# args[7] is the gtable -chr.file <- args[7] - -chr.size <- read.table(chr.file) - -# setting output files -r.output <- args[8] -overlap.output <- args[9] -npeaks.output <- args[10] -em.sav.output <- args[11] -uri.sav.output <- args[12] - -# sink(paste(output.prefix, "-Rout.txt", sep="")) -sink(r.output) - -############# process the data -cat("is.broadpeak", is.broadpeak, "\n") -# process data, summit: the representation of the location of summit -rep1 <- process.narrowpeak(paste(peakfile1, sep=""), chr.size, half.width=half.width, summit="offset", broadpeak=is.broadpeak) -rep2 <- process.narrowpeak(paste(peakfile2, sep=""), chr.size, half.width=half.width, summit="offset", broadpeak=is.broadpeak) - -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")) - -if(args[3]==-1){ - cat(paste("half.width=", "reported", "\n")) -}else{ - cat(paste("half.width=", half.width, "\n")) -} -cat(paste("significant measure=", sig.value, "\n")) - -# compute correspondence profile (URI) -uri.output <- compute.pair.uri(rep1$data.cleaned, rep2$data.cleaned, sig.value1=sig.value, sig.value2=sig.value, overlap.ratio=overlap.ratio) - -#uri.output <- compute.pair.uri(rep1$data.cleaned, rep2$data.cleaned) - -cat(paste("URI is done\n")) - -# save output -# save(uri.output, file=paste(output.prefix, "-uri.sav", sep="")) -save(uri.output, file=uri.sav.output) -cat(paste("URI is saved at: ", uri.sav.output)) - - -# EM procedure for inference -em.output <- fit.em(uri.output$data12.enrich, fix.rho2=T) - -#em.output <- fit.2copula.em(uri.output$data12.enrich, fix.rho2=T, "gaussian") - -cat(paste("EM is done\n\n")) - -save(em.output, file=em.sav.output) -cat(paste("EM is saved at: ", em.sav.output)) - - -# write em output into a file - -cat(paste("EM estimation for the following files\n", peakfile1, "\n", peakfile2, "\n", sep="")) - -print(em.output$em.fit$para) - -# add on 3-29-10 -# output both local idr and IDR -idr.local <- 1-em.output$em.fit$e.z -IDR <- c() -o <- order(idr.local) -IDR[o] <- cumsum(idr.local[o])/c(1:length(o)) - - -write.out.data <- data.frame(chr1=em.output$data.pruned$sample1[, "chr"], - start1=em.output$data.pruned$sample1[, "start.ori"], - stop1=em.output$data.pruned$sample1[, "stop.ori"], - sig.value1=em.output$data.pruned$sample1[, "sig.value"], - chr2=em.output$data.pruned$sample2[, "chr"], - start2=em.output$data.pruned$sample2[, "start.ori"], - stop2=em.output$data.pruned$sample2[, "stop.ori"], - sig.value2=em.output$data.pruned$sample2[, "sig.value"], - idr.local=1-em.output$em.fit$e.z, IDR=IDR) - -# write.table(write.out.data, file=paste(output.prefix, "-overlapped-peaks.txt", sep="")) -write.table(write.out.data, file=overlap.output) -cat(paste("Write overlapped peaks and local idr to: ", overlap.output, sep="")) - -# number of peaks passing IDR range (0.01-0.25) -IDR.cutoff <- seq(0.01, 0.25, by=0.01) -idr.o <- order(write.out.data$idr.local) -idr.ordered <- write.out.data$idr.local[idr.o] -IDR.sum <- cumsum(idr.ordered)/c(1:length(idr.ordered)) - -IDR.count <- c() -n.cutoff <- length(IDR.cutoff) -for(i in 1:n.cutoff){ - IDR.count[i] <- sum(IDR.sum <= IDR.cutoff[i]) -} - - -# write the number of peaks passing various IDR range into a file -idr.cut <- data.frame(peakfile1, peakfile2, IDR.cutoff=IDR.cutoff, IDR.count=IDR.count) -write.table(idr.cut, file=npeaks.output, append=T, quote=F, row.names=F, col.names=F) -cat(paste("Write number of peaks above IDR cutoff [0.01, 0.25]: ","npeaks-aboveIDR.txt\n", sep="")) - -mar.mean <- get.mar.mean(em.output$em.fit) - -cat(paste("Marginal mean of two components:\n")) -print(mar.mean) - -sink() - - |
| b |
| diff -r 48767bec000d -r 16bbe277d15d batch-consistency-plot.r --- a/batch-consistency-plot.r Thu Jan 17 15:45:48 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| @@ -1,68 +0,0 @@ -# 1-20-10 Qunhua Li -# -# This program first plots correspondence curve and IDR threshold plot -# (i.e. number of selected peaks vs IDR) for each pair of sample -# -# usage: -# Rscript batch-consistency-plot-merged.r [npairs] [output.dir] [input.file.prefix 1, 2, 3 ...] -# [npairs]: integer, number of consistency analyses -# (e.g. if 2 replicates, npairs=1, if 3 replicates, npairs=3 -# [output.prefix]: output directory and file name prefix for plot eg. /plots/idrPlot -# [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 -# - -args <- commandArgs(trailingOnly=T) -npair <- args[1] # number of curves to plot on the same figure -output.file.prefix <- args[2] # file name for plot, generated from script at the outer level -df.txt <- 10 -ntemp <- as.numeric(npair) -saved.file.prefix <- list() # identifier of filenames that contain the em and URI results -source("/mnt/galaxyTools/galaxy-central/tools/modENCODE_DCC_tools/idr/functions-all-clayton-12-13.r") - -uri.list <- list() -uri.list.match <- list() -ez.list <- list() -legend.txt <- c() -em.output.list <- list() -uri.output.list <- list() - -for(i in 1:npair){ - saved.file.prefix[i] <- args[2+i] - - load(paste(saved.file.prefix[i], "-uri.sav", sep="")) - load(paste(saved.file.prefix[i], "-em.sav", sep="")) - - uri.output.list[[i]] <- uri.output - em.output.list[[i]] <- em.output - - ez.list[[i]] <- get.ez.tt.all(em.output, uri.output.list[[i]]$data12.enrich$merge1, - uri.output.list[[i]]$data12.enrich$merge2) # reverse =T for error rate - - # URI for all peaks - uri.list[[i]] <- uri.output$uri.n - # URI for matched peaks - uri.match <- get.uri.matched(em.output$data.pruned, df=df.txt) - uri.list.match[[i]] <- uri.match$uri.n - - file.name <- unlist(strsplit(as.character(saved.file.prefix[i]), "/")) - - legend.txt[i] <- paste(i, "=", file.name[length(file.name)]) - -} - -plot.uri.file <- paste(output.file.prefix, "-plot.ps", sep="") - -############# plot and report output -# plot correspondence curve for each pair, -# plot number of selected peaks vs IDR -# plot all into 1 file -postscript(paste(output.file.prefix, "-plot.ps", sep="")) -par(mfcol=c(2,3), mar=c(5,6,4,2)+0.1) -plot.uri.group(uri.list, NULL, file.name=NULL, c(1:npair), title.txt="all peaks") -plot.uri.group(uri.list.match, NULL, file.name=NULL, c(1:npair), title.txt="matched peaks") -plot.ez.group(ez.list, plot.dir=NULL, file.name=NULL, legend.txt=c(1:npair), y.lim=c(0, 0.6)) -plot(0, 1, type="n", xlim=c(0,1), ylim=c(0,1), xlab="", ylab="", xaxt="n", yaxt="n") # legends -legend(0, 1, legend.txt, cex=0.6) - -dev.off() - |
| b |
| diff -r 48767bec000d -r 16bbe277d15d functions-all-clayton-12-13.r --- a/functions-all-clayton-12-13.r Thu Jan 17 15:45:48 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| b'@@ -1,3099 +0,0 @@\n-# revised on 2-20-10\n-# - fix error in pass.structure: reverse rank.combined, so that big sig.value\n-# are ranked with small numbers (1, 2, ...)\n-# - fix error on get.ez.tt.all: get ez.cutoff from sorted e.z\n-\n-#\n-# modified EM procedure to compute empirical CDF more precisely - 09/2009\n-\n-\n-\n-# this file contains the functions for \n-# 1. computing the correspondence profile (upper rank intersection and derivatives)\n-# 2. inference of copula mixture model\n-#\n-# It also has functions for\n-# 1. reading peak caller results\n-# 2. processing and matching called peaks\n-# 3. plotting results\n-\n-\n-################ read peak caller results\n-\n-# process narrow peak format\n-# some peak callers may not report q-values, p-values or fold of enrichment\n-# need further process before comparison\n-#\n-# stop.exclusive: Is the basepair of peak.list$stop exclusive? In narrowpeak and broadpeak format they are exclusive.\n-# If it is exclusive, we need subtract peak.list$stop by 1 to avoid the same basepair being both a start and a stop of two \n-# adjacent peaks, which creates trouble for finding correct intersect \n-process.narrowpeak <- function(narrow.file, chr.size, half.width=NULL, summit="offset", stop.exclusive=T, broadpeak=F){\n-\n-\n- aa <- read.table(narrow.file)\n-\n- if(broadpeak){\n- bb.ori <- data.frame(chr=aa$V1, start=aa$V2, stop=aa$V3, signal.value=aa$V7, p.value=aa$V8, q.value=aa$V9)\n- }else{\n- bb.ori <- data.frame(chr=aa$V1, start=aa$V2, stop=aa$V3, signal.value=aa$V7, p.value=aa$V8, q.value=aa$V9, summit=aa$V10)\n- }\n- \n- if(summit=="summit"){\n- bb.ori$summit <- bb.ori$summit-bb.ori$start # change summit to offset to avoid error when concatenating chromosomes\n- }\n- \n- bb <- concatenate.chr(bb.ori, chr.size)\n-\n- #bb <- bb.ori\n-\n- # remove the peaks that has the same start and stop value\n- bb <- bb[bb$start != bb$stop,]\n-\n- if(stop.exclusive==T){\n- bb$stop <- bb$stop-1\n- }\n-\n- if(!is.null(half.width)){\n- bb$start.ori <- bb$start \n- bb$stop.ori <- bb$stop \n-\n- # if peak is narrower than the specified window, stay with its width\n- # otherwise chop wider peaks to specified width\n- width <- bb$stop-bb$start +1\n- is.wider <- width > 2*half.width\n-\n- if(summit=="offset" | summit=="summit"){ # if summit is offset from start\n- bb$start[is.wider] <- bb$start.ori[is.wider] + bb$summit[is.wider]-half.width\n- bb$stop[is.wider] <- bb$start.ori[is.wider] + bb$summit[is.wider]+half.width\n- } else { \n- if(summit=="unknown"){\n- bb$start[is.wider] <- bb$start.ori[is.wider]+round(width[is.wider]/2) - half.width\n- bb$stop[is.wider] <- bb$start.ori[is.wider]+round(width[is.wider]/2) + half.width\n- }\n- }\n- }\n-\n- bb <- clean.data(bb)\n- invisible(list(data.ori=bb.ori, data.cleaned=bb))\n-}\n-\n-# clean data \n-# and concatenate chromosomes if needed\n-clean.data <- function(adata){\n-\n- # remove the peaks that has the same start and stop value\n- adata <- adata[adata$start != adata$stop,]\n-\n- # if some stops and starts are the same, need fix them\n- stop.in.start <- is.element(adata$stop, adata$start)\n- n.fix <- sum(stop.in.start)\n- if(n.fix >0){\n- print(paste("Fix", n.fix, "stops\\n"))\n- adata$stop[stop.in.start] <- adata$stop[stop.in.start]-1 \n- } \n- \n- return(adata) \n-}\n-\n-# concatenate peaks\n-# peaks: the dataframe to have all the peaks\n-# chr.file: the file to keep the length of each chromosome \n-# chr files should come from the species that the data is from\n-#concatenate.chr <- function(peaks, chr.size){\n-\n- # chr.size <- read.table(chr.file)\n-# chr.o <- order(chr.size[,1])\n-# chr.size <- chr.size[chr.o,]\n-#\n-# chr.shift <- cumsum(c(0, chr.size[-nrow(chr.size),2]))\n-# chr.size.cum <- data.frame(chr=chr.size[,1], shift=chr.shift) \n-#\n-# for(i in 1:nrow(chr.size)){\n-# is.in <- as.character(peaks$chr) == as.character(chr.size.cum$chr[i])\n-# if(sum(is.in)>0){\n-# peaks[is.in,]$start <- peaks[is.in,]$start + chr.size.cum$shi'..b'\n-# density of binomial distribution with equal mean and sigma on both dimensions\n-d.binormal <- function(z.1, z.2, mu, sigma, rho){\n-\n- loglik <- (-log(2)-log(pi)-2*log(sigma) - log(1-rho^2)/2 - (0.5/(1-rho^2)/sigma^2)*((z.1-mu)^2 -2*rho*(z.1-mu)*(z.2-mu) + (z.2-mu)^2))\n-\n- return(loglik)\n-}\n-\n-# E-step for computing the latent strucutre\n-# e.z is the prob to be in the consistent group\n-# e.step for estimating posterior prob\n-# z.1 and z.2 can be vectors or scalars\n-e.step.2normal <- function(z.1, z.2, mu, sigma, rho, p){\n-\n- e.z <- p/((1-p)*exp(d.binormal(z.1, z.2, 0, 1, 0)-d.binormal(z.1, z.2, mu, sigma, rho))+ p)\n- \n- invisible(e.z)\n-}\n-\n-# M-step for computing the latent structure\n-# m.step for estimating proportion, mean, sd and correlation coefficient\n-m.step.2normal <- function(z.1, z.2, e.z){\n-\n- p <- mean(e.z)\n- mu <- sum((z.1+z.2)*e.z)/2/sum(e.z) \n- sigma <- sqrt(sum(e.z*((z.1-mu)^2+(z.2-mu)^2))/2/sum(e.z))\n- rho <- 2*sum(e.z*(z.1-mu)*(z.2-mu))/(sum(e.z*((z.1-mu)^2+(z.2-mu)^2)))\n-\n- return(list(p=p, mu=mu, sigma=sigma, rho=rho))\n-}\n-\n-\n-# assume top p percent of observations are true\n-# x and y are ranks, estimate\n-init <- function(x, y, x.label){\n- \n- x.o <- order(x)\n-\n- x.ordered <- x[x.o]\n- y.ordered <- y[x.o]\n- x.label.ordered <- x.label[x.o]\n- \n- n <- length(x)\n- p <- sum(x.label)/n\n- \n- rho <- cor(x.ordered[1:ceiling(p*n)], y.ordered[1:ceiling(p*n)])\n-\n- temp <- find.mu.sigma(x.ordered, x.label.ordered)\n- mu <- temp$mu\n- sigma <- temp$sigma\n- \n- invisible(list(mu=mu, sigma=sigma, rho=rho, p=p))\n-\n-}\n-\n-# find mu and sigma if the distributions of marginal ranks are known\n-# take the medians of the two dist and map back to the original\n-init.dist <- function(f0, f1){\n-\n- # take the median in f0\n- index.median.0 <- which(f0$cdf>0.5)[1]\n- q.0.small <- f0$cdf[index.median.0] # because f0 and f1 have the same bins\n- q.1.small <- f1$cdf[index.median.0]\n-\n- # take the median in f1\n- index.median.1 <- which(f1$cdf>0.5)[1]\n- q.0.big <- f0$cdf[index.median.1] # because f0 and f1 have the same bins\n- q.1.big <- f1$cdf[index.median.1]\n-\n- # find pseudo value for x.middle[1] on normal(0,1) \n- pseudo.small.0 <- qnorm(q.0.small, mean=0, sd=1)\n- pseudo.small.1 <- qnorm(q.1.small, mean=0, sd=1)\n-\n- # find pseudo value for x.middle[2] on normal(0,1) \n- pseudo.big.0 <- qnorm(q.0.big, mean=0, sd=1)\n- pseudo.big.1 <- qnorm(q.1.big, mean=0, sd=1)\n-\n- mu <- (pseudo.small.0*pseudo.big.1 - pseudo.small.1*pseudo.big.0)/(pseudo.big.1-pseudo.small.1) \n-\n- sigma <- (pseudo.small.0-mu)/pseudo.small.1\n-\n- return(list(mu=mu, sigma=sigma)) \n-}\n-\n-# generate labels\n-\n-# find the part of data with overlap\n-\n-# find the percentile on noise and signal\n-\n-# Suppose there are signal and noise components, with mean=0 and sd=1 for noise\n-# x and x.label are the rank of the observations and their labels,\n-# find the mean and sd of the other component\n-# x.label takes values of 0 and 1\n-find.mu.sigma <- function(x, x.label){\n-\n- x.0 <- x[x.label==0]\n- x.1 <- x[x.label==1]\n-\n- n.x0 <- length(x.0)\n- n.x1 <- length(x.1)\n-\n- x.end <- c(min(x.0), min(x.1), max(x.0), max(x.1))\n- o <- order(x.end)\n- x.middle <- x.end[o][c(2,3)]\n-\n- # the smaller end of the overlap\n- q.1.small <- mean(x.1 <= x.middle[1])*n.x1/(n.x1+1)\n- q.0.small <- mean(x.0 <= x.middle[1])*n.x0/(n.x0+1)\n-\n- # the bigger end of the overlap\n- q.1.big <- mean(x.1 <= x.middle[2])*n.x1/(n.x1+1)\n- q.0.big <- mean(x.0 <= x.middle[2])*n.x0/(n.x0+1)\n-\n- # find pseudo value for x.middle[1] on normal(0,1) \n- pseudo.small.0 <- qnorm(q.0.small, mean=0, sd=1)\n- pseudo.small.1 <- qnorm(q.1.small, mean=0, sd=1)\n-\n- # find pseudo value for x.middle[2] on normal(0,1) \n- pseudo.big.0 <- qnorm(q.0.big, mean=0, sd=1)\n- pseudo.big.1 <- qnorm(q.1.big, mean=0, sd=1)\n-\n- mu <- (pseudo.small.0*pseudo.big.1 - pseudo.small.1*pseudo.big.0)/(pseudo.big.1-pseudo.small.1) \n-\n- sigma <- (pseudo.small.0-mu)/pseudo.small.1\n-\n- return(list(mu=mu, sigma=sigma))\n-}\n' |
| b |
| diff -r 48767bec000d -r 16bbe277d15d genome_tables/genome_table.dmel.r5.32.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/genome_tables/genome_table.dmel.r5.32.txt Thu Jan 17 15:57:35 2013 -0500 |
| b |
| @@ -0,0 +1,15 @@ +YHet 347038 +dmel_mitochondrion_genome 19517 +2L 23011544 +X 22422827 +3L 24543557 +4 1351857 +2R 21146708 +3R 27905053 +Uextra 29004656 +2RHet 3288761 +2LHet 368872 +3LHet 2555491 +3RHet 2517507 +U 10049037 +XHet 204112 |
| b |
| diff -r 48767bec000d -r 16bbe277d15d genome_tables/genome_table.human.hg18.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/genome_tables/genome_table.human.hg18.txt Thu Jan 17 15:57:35 2013 -0500 |
| b |
| @@ -0,0 +1,25 @@ +chr1 247249719 +chr2 242951149 +chr3 199501827 +chr4 191273063 +chr5 180857866 +chr6 170899992 +chr7 158821424 +chr8 146274826 +chr9 140273252 +chr10 135374737 +chr11 134452384 +chr12 132349534 +chr13 114142980 +chr14 106368585 +chr15 100338915 +chr16 88827254 +chr17 78774742 +chr18 76117153 +chr19 63811651 +chr20 62435964 +chr21 46944323 +chr22 49691432 +chrX 154913754 +chrY 57772954 +chrM 16571 |
| b |
| diff -r 48767bec000d -r 16bbe277d15d genome_tables/genome_table.human.hg19.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/genome_tables/genome_table.human.hg19.txt Thu Jan 17 15:57:35 2013 -0500 |
| b |
| @@ -0,0 +1,25 @@ +chr1 249250621 +chr2 243199373 +chr3 198022430 +chr4 191154276 +chr5 180915260 +chr6 171115067 +chr7 159138663 +chr8 146364022 +chr9 141213431 +chr10 135534747 +chr11 135006516 +chr12 133851895 +chr13 115169878 +chr14 107349540 +chr15 102531392 +chr16 90354753 +chr17 81195210 +chr18 78077248 +chr19 59128983 +chr20 63025520 +chr21 48129895 +chr22 51304566 +chrX 155270560 +chrY 59373566 +chrM 16571 |
| b |
| diff -r 48767bec000d -r 16bbe277d15d genome_tables/genome_table.mm9.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/genome_tables/genome_table.mm9.txt Thu Jan 17 15:57:35 2013 -0500 |
| b |
| @@ -0,0 +1,22 @@ +chr1 197195432 +chr2 181748087 +chr3 159599783 +chr4 155630120 +chr5 152537259 +chr6 149517037 +chr7 152524553 +chr8 131738871 +chr9 124076172 +chr10 129993255 +chr11 121843856 +chr12 121257530 +chr13 120284312 +chr14 125194864 +chr15 103494974 +chr16 98319150 +chr17 95272651 +chr18 90772031 +chr19 61342430 +chrX 166650296 +chrY 15902555 +chrM 16299 |
| b |
| diff -r 48767bec000d -r 16bbe277d15d genome_tables/genome_table.worm.ws220.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/genome_tables/genome_table.worm.ws220.txt Thu Jan 17 15:57:35 2013 -0500 |
| b |
| @@ -0,0 +1,7 @@ +I 15072423 +II 15279345 +III 13783700 +IV 17493793 +V 20924149 +X 17718866 +MtDNA 13794 \ No newline at end of file |
| b |
| diff -r 48767bec000d -r 16bbe277d15d idrPlotDef.xml --- a/idrPlotDef.xml Thu Jan 17 15:45:48 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| @@ -1,31 +0,0 @@ -<!-- -Usage of THIS SCRIPT: ./idrPlotWrapper.sh em uri outputfile ---> - -<tool id="batch_consistency_plot" name="IDR-Plot"> - <description>Plot Consistency Analysis on IDR output files</description> - <command interpreter="bash">idrPlotWrapper.sh $em $uri $outputfile</command> - <inputs> - <param name="uri" type="data" label="uri.sav file (output from IDR consistency analysis)"/> - <param name="em" type="data" label="em.sav file (output from IDR consistency analysis)"/> - </inputs> - <outputs> - <data format="pdf" name="outputfile" label="IDR-plot.pdf"/> - </outputs> - - <tests> - <test> -<!-- - <param name="input" value="fa_gc_content_input.fa"/> - <output name="out_file1" file="fa_gc_content_output.txt"/> ---> - </test> - </tests> - - <help> -Plots the correspondence curve and IDR threshold (i.e. number of selected peaks vs IDR) for each pair of samples. -Uses the output of IDR consistency analysis and produces a downloadable PDF file containing the graphical analysis. -This is a part of the IDR package. For more information about IDR, see https://sites.google.com/site/anshulkundaje/projects/idr. - </help> - -</tool> |
| b |
| diff -r 48767bec000d -r 16bbe277d15d idrPlotWrapper.sh --- a/idrPlotWrapper.sh Thu Jan 17 15:45:48 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| @@ -1,39 +0,0 @@ -#!/bin/bash - -# idrPlotWrapper.sh -# OICR: Kar Ming Chu -# July 2012 - -# BASH wrapper for batch-consistency-plot.r (part of the IDR package) -# For use with Galaxy - -# Usage of batch-consistency-plot.r: Rscript batch-consistency-plot-merged.r [npairs] [output.dir] [input.file.prefix 1, 2, 3 ...] -# npairs - will be a constant, since Galaxy requires explicit control over input and output files - -# Usage of THIS SCRIPT: ./idrPlotWrapper.sh em uri outputfile -# em - em.sav file provided by Galaxy -# uri - uri.sav file provided by Galaxy -# outputfile - output file name specified by Galaxy - -main() { - EM="${1}" # absolute file path to em.sav file, provided by Galaxy - URI="${2}" # absolute file parth to uri.sav file, provided by Galaxy - OUTFILE="${3}" # name of desired output file - - cp "${EM}" ./idrPlot-em.sav # cp to this directory and rename so they can be found by idrPlot - cp "${URI}" ./idrPlot-uri.sav - - Rscript /mnt/galaxyTools/galaxy-central/tools/modENCODE_DCC_tools/idr/batch-consistency-plot.r 1 ./idrPlot idrPlot - - # convert post script to pdf file - ps2pdf ./idrPlot-plot.ps ./idrPlot-plot.pdf - - # rename to output file name - mv ./idrPlot-plot.pdf "${OUTFILE}" - - # clean up - rm idrPlot-em.sav - rm idrPlot-uri.sav -} - -main "${@}" |
| b |
| diff -r 48767bec000d -r 16bbe277d15d idrToolDef.xml --- a/idrToolDef.xml Thu Jan 17 15:45:48 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| @@ -1,61 +0,0 @@ -<!-- - -Script Usage: -Rscript batch-consistency-analysis.r -../3066_rep1_VS_input0.macs14.out.regionPeak -../3066_rep2_VS_input0.macs14.out.regionPeak -1000 -3066_rep1_VS_rep2 -0 -F -p.value -genome_table.txt [ drop down to select ] ---> - -<tool id="batch_consistency_analysis_2" name="IDR"> - <description>Consistency Analysis on a pair of narrowPeak files</description> - <command interpreter="Rscript">batch-consistency-analysis.r $input1 $input2 $halfwidth $overlap $option $sigvalue $gtable $rout $aboveIDR $ratio $emSav $uriSav</command> - <inputs> - <param format="narrowPeak" name="input1" type="data" label="First NarrowPeak File"/> - <param format="narrowPeak" name="input2" type="data" label="Second NarrowPeak File"/> - <param name="halfwidth" size="4" type="integer" value="1000" label="Half-Width" help="-1 if using reported peak width"/> -<!-- <param name="outputprefix" type="text" size="50" label="Output Prefix" value="3066_rep1_VS_rep2"/> --> - <param name="option" type="select" label="File Type" value="F"> - <option value="F">Narrow Peak</option> - <option value="T">Broad Peak</option> - </param> - <param name="overlap" size="4" type="float" value="0" label="Over-Lap Ratio" help="Between 0 and 1, inclusively" min="0" max="1"/> - <param name="sigvalue" type="select" label="Significant Value" value="p.value" help="Select p-value if the input peak files are generated by MAC. Select q-value if the input peak files are generated by SPP."> - <option value="p.value">p-value</option> - <option value="q.value">q-value</option> - <option value="signal.value">Significant Value</option> - </param> - <param name="gtable" type="select" label="Genome Table" value="/mnt/galaxyTools/galaxy-central/tools/modENCODE_DCC_tools/idr/genome_tables/genome_table.worm.ws220.txt"> - <option value="/mnt/galaxyTools/galaxy-central/tools/modENCODE_DCC_tools/idr/genome_tables/genome_table.human.hg19.txt">human hg19</option> - <option value="/mnt/galaxyTools/galaxy-central/tools/modENCODE_DCC_tools/idr/genome_tables/genome_table.mm9.txt">mouse mm9</option> - <option value="/mnt/galaxyTools/galaxy-central/tools/modENCODE_DCC_tools/idr/genome_tables/genome_table.worm.ws220.txt">worm ws220</option> - <option value="/mnt/galaxyTools/galaxy-central/tools/modENCODE_DCC_tools/idr/genome_tables/genome_table.dmel.r5.32.txt">dmel r5.32</option> - </param> - </inputs> - <outputs> - <data format="txt" name="rout" label="IDR.Rout.txt"/> - <data format="txt" name="aboveIDR" label="IDR.npeaks-aboveIDR.txt"/> - <data format="txt" name="ratio" label="IDR.overlapped-peaks.txt"/> - <data format="txt" name="emSav" label="IDR.em.sav"/> - <data format="txt" name="uriSav" label="IDR.uri.sav"/> - </outputs> - - <tests> - <test> -<!-- - <param name="input" value="fa_gc_content_input.fa"/> - <output name="out_file1" file="fa_gc_content_output.txt"/> ---> - </test> - </tests> - - <help> -Reproducibility is essential to reliable scientific discovery in high-throughput experiments. The IDR (Irreproducible Discovery Rate) framework is a unified approach to measure the reproducibility of findings identified from replicate experiments and provide highly stable thresholds based on reproducibility. Unlike the usual scalar measures of reproducibility, the IDR approach creates a curve, which quantitatively assesses when the findings are no longer consistent across replicates. In layman's terms, the IDR method compares a pair of ranked lists of identifications (such as ChIP-seq peaks). These ranked lists should not be pre-thresholded i.e. they should provide identifications across the entire spectrum of high confidence/enrichment (signal) and low confidence/enrichment (noise). The IDR method then fits the bivariate rank distributions over the replicates in order to separate signal from noise based on a defined confidence of rank consistency and reproducibility of identifications i.e the IDR threshold. For more information on IDR, see https://sites.google.com/site/anshulkundaje/projects/idr - </help> - -</tool> |