| Previous changeset 12:8781805e48e7 (2013-01-17) Next changeset 14:58a1d7ef98ef (2013-01-21) |
|
Commit message:
Uploaded |
|
added:
genome_tables/.DS_Store genome_tables/._.DS_Store |
|
removed:
batch-consistency-analysis.r batch-consistency-plot.r functions-all-clayton-12-13.r idrPlotDef.xml idrPlotWrapper.sh idrToolDef.xml |
| b |
| diff -r 8781805e48e7 -r 2ebe483e9895 batch-consistency-analysis.r --- a/batch-consistency-analysis.r Thu Jan 17 15:59:30 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 8781805e48e7 -r 2ebe483e9895 batch-consistency-plot.r --- a/batch-consistency-plot.r Thu Jan 17 15:59:30 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 8781805e48e7 -r 2ebe483e9895 functions-all-clayton-12-13.r --- a/functions-all-clayton-12-13.r Thu Jan 17 15:59:30 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 8781805e48e7 -r 2ebe483e9895 genome_tables/.DS_Store |
| b |
| Binary file genome_tables/.DS_Store has changed |
| b |
| diff -r 8781805e48e7 -r 2ebe483e9895 genome_tables/._.DS_Store |
| b |
| Binary file genome_tables/._.DS_Store has changed |
| b |
| diff -r 8781805e48e7 -r 2ebe483e9895 idrPlotDef.xml --- a/idrPlotDef.xml Thu Jan 17 15:59:30 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 8781805e48e7 -r 2ebe483e9895 idrPlotWrapper.sh --- a/idrPlotWrapper.sh Thu Jan 17 15:59:30 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 8781805e48e7 -r 2ebe483e9895 idrToolDef.xml --- a/idrToolDef.xml Thu Jan 17 15:59:30 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> |