Mercurial > repos > ziru-zhou > spp_phantompeak

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/run_spp.R	Tue Dec 18 10:44:50 2012 -0500
@@ -0,0 +1,885 @@
+# run_spp.R
+# =============
+# Author: Anshul Kundaje, Computer Science Dept., Stanford University
+# Email: akundaje@stanford.edu
+# Last updated: Oct 8, 2010
+# =============
+# MANDATORY ARGUMENTS
+# -c=<ChIP_tagAlign/BAMFile>, full path and name of tagAlign/BAM file (can be gzipped) (FILE EXTENSION MUST BE tagAlign.gz, tagAlign, bam or bam.gz)
+# MANDATORY ARGUMENT FOR PEAK CALLING
+# -i=<Input_tagAlign/BAMFile>, full path and name of tagAlign/BAM file (can be gzipped) (FILE EXTENSION MUST BE tagAlign.gz, tagAlign, bam or bam.gz)
+# OPTIONAL ARGUMENTS
+# -s=<min>:<step>:<max> , strand shifts at which cross-correlation is evaluated, default=-100:5:600
+# -speak=<strPeak>, user-defined cross-correlation peak strandshift
+# -x=<min>:<max>, strand shifts to exclude (This is mainly to avoid phantom peaks) default=10:(readlen+10)
+# -p=<nodes> , number of parallel processing nodes, default=NULL
+# -fdr=<falseDisoveryRate> , false discovery rate threshold for peak calling
+# -npeak=<numPeaks>, threshold on number of peaks to call
+# -tmpdir=<tempdir> , Temporary directory (if not specified R function tempdir() is used)
+# -filtchr=<chrnamePattern> , Pattern to use to remove tags that map to specific chromosomes e.g. _ will remove all tags that map to chromosomes with _ in their name
+# OUTPUT PARAMETERS
+# -odir=<outputDirectory> name of output directory (If not set same as ChIP file directory is used)
+# -savn=<narrowpeakfilename> OR -savn NarrowPeak file name
+# -savr=<regionpeakfilename> OR -savr RegionPeak file name
+# -savd=<rdatafile> OR -savd , save Rdata file
+# -savp=<plotdatafile> OR -savp , save cross-correlation plot
+# -out=<resultfile>, append peakshift result to a file
+#      format:Filename<tab>numReads<tab>estFragLen<tab>corr_estFragLen<tab>PhantomPeak<tab>corr_phantomPeak<tab>argmin_corr<tab>min_corr<tab>phantomPeakCoef<tab>relPhantomPeakCoef<tab>QualityTag
+# -rf , if plot or rdata or narrowPeak file exists replace it. If not used then the run is aborted if the plot or Rdata or narrowPeak file exists
+# -clean, if present will remove the original chip and control files after reading them in. CAUTION: Use only if the script calling run_spp.R is creating temporary files
+
+args <- commandArgs(trailingOnly=TRUE); # Read Arguments from command line
+nargs = length(args); # number of arguments
+
+# ###########################################################################
+# AUXILIARY FUNCTIONS
+# ###########################################################################
+
+print.usage <- function() {
+# ===================================
+# Function will print function usage
+# ===================================
+	cat('Usage: Rscript run_spp.R <options>\n',file=stderr())
+	cat('MANDATORY ARGUMENTS\n',file=stderr())
+	cat('-c=<ChIP_alignFile>, full path and name (or URL) of tagAlign/BAM file (can be gzipped)(FILE EXTENSION MUST BE tagAlign.gz, tagAlign, bam or bam.gz) \n',file=stderr())
+	cat('MANDATORY ARGUMENTS FOR PEAK CALLING\n',file=stderr())
+	cat('-i=<Input_alignFile>, full path and name (or URL) of tagAlign/BAM file (can be gzipped) (FILE EXTENSION MUST BE tagAlign.gz, tagAlign, bam or bam.gz) \n',file=stderr())
+	cat('OPTIONAL ARGUMENTS\n',file=stderr())
+	cat('-s=<min>:<step>:<max> , strand shifts at which cross-correlation is evaluated, default=-100:5:600\n',file=stderr())
+	cat('-speak=<strPeak>, user-defined cross-correlation peak strandshift\n',file=stderr())
+	cat('-x=<min>:<max>, strand shifts to exclude (This is mainly to avoid region around phantom peak) default=10:(readlen+10)\n',file=stderr())
+	cat('-p=<nodes> , number of parallel processing nodes, default=0\n',file=stderr())
+	cat('-fdr=<falseDisoveryRate> , false discovery rate threshold for peak calling\n',file=stderr())
+	cat('-npeak=<numPeaks>, threshold on number of peaks to call\n',file=stderr())
+	cat('-tmpdir=<tempdir> , Temporary directory (if not specified R function tempdir() is used)\n',file=stderr())
+	cat('-filtchr=<chrnamePattern> , Pattern to use to remove tags that map to specific chromosomes e.g. _ will remove all tags that map to chromosomes with _ in their name\n',file=stderr())
+	cat('OUTPUT ARGUMENTS\n',file=stderr())
+	cat('-odir=<outputDirectory> name of output directory (If not set same as ChIP file directory is used)\n',file=stderr())
+	cat('-savn=<narrowpeakfilename> OR -savn NarrowPeak file name (fixed width peaks)\n',file=stderr())
+	cat('-savr=<regionpeakfilename> OR -savr RegionPeak file name (variable width peaks with regions of enrichment)\n',file=stderr())
+	cat('-savd=<rdatafile> OR -savd, save Rdata file\n',file=stderr())
+	cat('-savp=<plotdatafile> OR -savp, save cross-correlation plot\n',file=stderr())
+	cat('-out=<resultfile>, append peakshift/phantomPeak results to a file\n',file=stderr())
+	cat('     format:Filename<tab>numReads<tab>estFragLen<tab>corr_estFragLen<tab>PhantomPeak<tab>corr_phantomPeak<tab>argmin_corr<tab>min_corr<tab>phantomPeakCoef<tab>relPhantomPeakCoef<tab>QualityTag)\n',file=stderr())
+	cat('-rf, if plot or rdata or narrowPeak file exists replace it. If not used then the run is aborted if the plot or Rdata or narrowPeak file exists\n',file=stderr())
+	cat('-clean, if present will remove the original chip and control files after reading them in. CAUTION: Use only if the script calling run_spp.R is creating temporary files\n',file=stderr())
+} # end: print.usage()
+
+get.file.parts <- function(file.fullpath) {
+# ===================================
+# Function will take a file name with path and split the file name into
+# path, fullname, name and ext
+# ===================================
+	if (! is.character(file.fullpath)) {
+		stop('File name must be a string')
+	}
+
+	file.parts <- strsplit(as.character(file.fullpath), .Platform$file.sep, fixed=TRUE)[[1]] # split on file separator
+
+	if (length(file.parts) == 0) { # if empty file name
+		return(list(path='',
+						fullname='',
+						name='',
+						ext='')
+		)
+	} else {
+		if (length(file.parts) == 1) { # if no path then just the file name itself
+			file.path <- '.'
+			file.fullname <- file.parts
+		} else {
+			file.path <- paste(file.parts[1:(length(file.parts)-1)], collapse=.Platform$file.sep) # 1:last-1 token is path
+			file.fullname <- file.parts[length(file.parts)] # last token is filename
+		}
+		file.fullname.parts <- strsplit(file.fullname,'.',fixed=TRUE)[[1]] # split on .
+		if (length(file.fullname.parts) == 1) { # if no extension
+			file.ext <- ''
+			file.name <- file.fullname.parts
+		} else {
+			file.ext <- paste('.', file.fullname.parts[length(file.fullname.parts)], sep="") # add the . to the last token
+			file.name <- paste(file.fullname.parts[1:(length(file.fullname.parts)-1)], collapse=".")
+		}
+		return(list(path=file.path,
+						fullname=file.fullname,
+						name=file.name,
+						ext=file.ext))
+	}
+} # end: get.file.parts()
+
+parse.arguments <- function(args) {
+# ===================================
+# Function will parse arguments
+# ===================================
+	# Set arguments to default values
+	chip.file <- NA  # main ChIP tagAlign/BAM file name
+	isurl.chip.file <- FALSE # flag indicating whether ChIP file is a URL
+	control.file <- NA # control tagAlign/BAM file name
+	isurl.control.file <- FALSE # flag indicating whether control file is a URL
+	sep.min <- -100  # min strand shift
+	sep.max <- 600  # max strand shift
+	sep.bin <- 5    # increment for strand shift
+	sep.peak <- NA # user-defined peak shift
+	exclude.min <- 10 # lowerbound of strand shift exclusion region
+	exclude.max <- NaN # upperbound of strand shift exclusion region
+	n.nodes <- NA # number of parallel processing nodes
+	fdr <- 0.01 # false discovery rate threshold for peak calling
+	npeak <- NA # threshold on number of peaks to call
+	temp.dir <- tempdir() # temporary directory
+	chrname.rm.pattern <- NA # chromosome name pattern used to remove tags
+	output.odir <- NA # Output directory name
+	output.npeak.file <- NA # Output narrowPeak file name
+	output.rpeak.file <- NA # Output regionPeak file name
+	output.rdata.file <- NA # Rdata file
+	output.plot.file <- NA  # cross correlation plot file
+	output.result.file <- NA # result file
+	replace.flag <- FALSE # replace file flag
+	clean.files.flag <- FALSE # file deletion flag
+
+	# Parse arguments
+	for (each.arg in args) {
+
+		if (grepl('^-c=',each.arg)) { #-c=<chip.file>
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				chip.file <- arg.split[2] # second part is chip.file
+			} else {
+				stop('No tagAlign/BAM file name provided for parameter -c=')
+			}
+
+		} else if (grepl('^-i=',each.arg)) { #-i=<control.file>
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				control.file <- arg.split[2] # second part is control.file
+			} else {
+				stop('No tagAlign/BAM file name provided for parameter -i=')
+			}
+
+		} else if (grepl('^-s=',each.arg)) { #-s=<sep.min>:<sep.bin>:<sep.max>
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				sep.vals <- arg.split[2] # second part is sepmin:sepbin:sepmax
+				sep.vals.split <- strsplit(sep.vals,':',fixed=TRUE)[[1]] # split on :
+				if (length(sep.vals.split) != 3) { # must have 3 parts
+					stop('Strand shift limits must be specified as -s=sepmin:sepbin:sepmax')
+				} else {
+					if (any(is.na(as.numeric(sep.vals.split)))) { # check that sep vals are numeric
+						stop('Strand shift limits must be numeric values')
+					}
+					sep.min <- round(as.numeric(sep.vals.split[1]))
+					sep.bin <- round(as.numeric(sep.vals.split[2]))
+					sep.max <- round(as.numeric(sep.vals.split[3]))
+					if ((sep.min > sep.max) || (sep.bin > (sep.max - sep.min)) || (sep.bin < 0)) {
+						stop('Illegal separation values -s=sepmin:sepbin:sepmax')
+					}
+				}
+			} else {
+				stop('Strand shift limits must be specified as -s=sepmin:sepbin:sepmax')
+			}
+
+		} else if (grepl('^-speak=',each.arg)) { #-speak=<sep.peak> , user-defined cross-correlation peak strandshift
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				sep.peak <- arg.split[2] # second part is <sep.peak>
+				if (is.na(as.numeric(sep.peak))) { # check that sep.peak is numeric
+					stop('-speak=<sep.peak>: User defined peak shift must be numeric')
+				}
+				sep.peak <- as.numeric(sep.peak)
+			} else {
+				stop('User defined peak shift must be provided as -speak=<sep.peak>')
+			}
+
+		} else if (grepl('^-x=',each.arg)) { #-x=<exclude.min>:<exclude.max>
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				exclude.vals <- arg.split[2] # second part is excludemin:excludemax
+				exclude.vals.split <- strsplit(exclude.vals,':',fixed=TRUE)[[1]] # split on :
+				if (length(exclude.vals.split) != 2) { # must have 2 parts
+					stop('Exclusion limits must be specified as -x=excludemin:excludemax')
+				} else {
+					if (any(is.na(as.numeric(exclude.vals.split)))) { # check that exclude vals are numeric
+						stop('Exclusion limits must be numeric values')
+					}
+					exclude.min <- round(as.numeric(exclude.vals.split[1]))
+					exclude.max <- round(as.numeric(exclude.vals.split[2]))
+					if (exclude.min > exclude.max) {
+						stop('Illegal exclusion limits -x=excludemin:excludemax')
+					}
+				}
+			} else {
+				stop('Exclusion limits must be specified as -x=excludemin:excludemax')
+			}
+
+		} else if (grepl('^-p=',each.arg)) { #-p=<n.nodes> , number of parallel processing nodes, default=NULL
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				n.nodes <- arg.split[2] # second part is numnodes
+				if (is.na(as.numeric(n.nodes))) { # check that n.nodes is numeric
+					stop('-p=<numnodes>: numnodes must be numeric')
+				}
+				n.nodes <- round(as.numeric(n.nodes))
+			} else {
+				stop('Number of parallel nodes must be provided as -p=<numnodes>')
+			}
+
+		} else if (grepl('^-fdr=',each.arg)) { #-fdr=<fdr> , false discovery rate, default=0.01
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				fdr <- arg.split[2] # second part is fdr
+				if (is.na(as.numeric(fdr))) { # check that fdr is numeric
+					stop('-fdr=<falseDiscoveryRate>: false discovery rate must be numeric')
+				}
+				fdr <- as.numeric(fdr)
+			} else {
+				stop('False discovery rate must be provided as -fdr=<fdr>')
+			}
+
+		} else if (grepl('^-npeak=',each.arg)) { #-npeak=<numPeaks> , number of peaks threshold, default=NA
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				npeak <- arg.split[2] # second part is npeak
+				if (is.na(as.numeric(npeak))) { # check that npeak is numeric
+					stop('-npeak=<numPeaks>: threshold on number of peaks must be numeric')
+				}
+				npeak <- round(as.numeric(npeak))
+			} else {
+				stop('Threshold on number of peaks must be provided as -npeak=<numPeaks>')
+			}
+
+		} else if (grepl('^-tmpdir=',each.arg)) { #-tmpdir=<temp.dir>
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				temp.dir <- arg.split[2] # second part is temp.dir
+			} else {
+				stop('No temporary directory provided for parameter -tmpdir=')
+			}
+
+		} else if (grepl('^-filtchr=',each.arg)) { #-filtchr=<chrname.rm.pattern>
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				chrname.rm.pattern <- arg.split[2] # second part is chrname.rm.pattern
+			} else {
+				stop('No pattern provided for parameter -filtchr=')
+			}
+
+		} else if (grepl('^-odir=',each.arg)) { #-odir=<output.odir>
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				output.odir <- arg.split[2] # second part is output.odir
+			} else {
+				stop('No output directory provided for parameter -odir=')
+			}
+
+		} else if (grepl('^-savn',each.arg)) { # -savn=<output.npeak.file> OR -savn , save narrowpeak
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2])) {
+				output.npeak.file <- arg.split[2] #-savn=
+			} else if (each.arg=='-savn') {
+				output.npeak.file <- NULL # NULL indicates get the name from the main file name
+			} else {
+				stop('Argument for saving narrowPeak file must be -savn or -savn=<filename>')
+			}
+
+		} else if (grepl('^-savr',each.arg)) { # -savr=<output.rpeak.file> OR -savr , save regionpeak
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2])) {
+				output.rpeak.file <- arg.split[2] #-savr=
+			} else if (each.arg=='-savr') {
+				output.rpeak.file <- NULL # NULL indicates get the name from the main file name
+			} else {
+				stop('Argument for saving regionPeak file must be -savr or -savr=<filename>')
+			}
+
+		} else if (grepl('^-savd',each.arg)) { # -savd=<output.rdata.file> OR -savd , save Rdata file
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2])) {
+				output.rdata.file <- arg.split[2] #-savd=
+			} else if (each.arg=='-savd') {
+				output.rdata.file <- NULL # NULL indicates get the name from the main file name
+			} else {
+				stop('Argument for saving Rdata file must be -savd or -savd=<filename>')
+			}
+
+		} else if (grepl('^-savp',each.arg)) { # -savp=<output.plot.file> OR -savp , save cross-correlation plot
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2])) {
+				output.plot.file <- arg.split[2] #-savp=
+			} else if (each.arg=='-savp') {
+				output.plot.file <- NULL # NULL indicates get the name from the main file name
+			} else {
+				stop('Argument for saving Rdata file must be -savp or -savp=<filename>')
+			}
+
+		} else if (grepl('^-out=',each.arg)) { #-out=<output.result.file>
+
+			arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+			if (! is.na(arg.split[2]) ) {
+				output.result.file <- arg.split[2] # second part is output.result.file
+			} else {
+				stop('No result file provided for parameter -out=')
+			}
+
+		} else if (each.arg == '-rf') {
+
+			replace.flag <- TRUE
+
+		} else if (each.arg == '-clean') {
+
+			clean.files.flag <- TRUE
+
+		} else {
+
+			stop('Illegal argument ',each.arg)
+		}
+	}
+	# End: for loop
+
+	# Check mandatory arguments
+	if (is.na(chip.file)) {
+		stop('-c=<tagAlign/BAMFileName> is a mandatory argument')
+	}
+
+	if (is.na(control.file) && ! is.na(output.npeak.file)) {
+		stop('-i=<tagAlign/BAMFileName> is required for peak calling')
+	}
+
+	# Check if ChIP and control files are URLs
+	if (grepl('^http://',chip.file)) {
+		isurl.chip.file <- TRUE
+	}
+	if (grepl('^http://',control.file)) {
+		isurl.control.file <- TRUE
+	}
+
+	# If ChIP file is a URL output.odir MUST be specified
+	if (isurl.chip.file && is.na(output.odir)) {
+		stop('If ChIP file is a URL, then output directory MUST be specified')
+	}
+
+	# Check that ChIP and control files exist
+	if (isurl.chip.file) {
+		if (system(paste('wget -q --spider',chip.file)) != 0) {
+			stop('ChIP file URL not valid: ',chip.file)
+		}
+	} else if (!file.exists(chip.file)) {
+		stop('ChIP File:',chip.file,' does not exist')
+	}
+
+	if (!is.na(control.file)) {
+		if (isurl.control.file) {
+			if (system(paste('wget -q --spider',control.file)) != 0) {
+				stop('Control file URL not valid: ',control.file)
+			}
+		} else if (!file.exists(control.file)) {
+			stop('Control File:',control.file,' does not exist')
+		}
+	}
+
+	# Correct other arguments
+	if (is.na(output.odir)) { # Reconstruct output.odir if not provided
+		output.odir <- get.file.parts(chip.file)$path
+	}
+
+	if (is.null(output.npeak.file)) { # Reconstruct output.npeak.file if NULL
+		output.npeak.file <- file.path(output.odir, paste(get.file.parts(chip.file)$name, '_VS_', get.file.parts(control.file)$name,'.narrowPeak', sep=""))
+	}
+
+	if (is.null(output.rpeak.file)) { # Reconstruct output.rpeak.file if NULL
+		output.rpeak.file <- file.path(output.odir, paste(get.file.parts(chip.file)$name, '_VS_', get.file.parts(control.file)$name,'.regionPeak', sep=""))
+	}
+
+	if (is.null(output.rdata.file)) { # Reconstruct output.rdata.file if NULL
+		output.rdata.file <- file.path(output.odir, paste(get.file.parts(chip.file)$name, '.Rdata', sep=""))
+	}
+
+	if (is.null(output.plot.file)) { # Reconstruct output.plot.file if NULL
+		output.plot.file <- file.path(output.odir, paste(get.file.parts(chip.file)$name, '.pdf', sep=""))
+	}
+
+	return(list(chip.file=chip.file,
+					isurl.chip.file=isurl.chip.file,
+					control.file=control.file,
+					isurl.control.file=isurl.control.file,
+					sep.range=c(sep.min,sep.bin,sep.max),
+					sep.peak=sep.peak,
+					ex.range=c(exclude.min,exclude.max),
+					n.nodes=n.nodes,
+					fdr=fdr,
+					npeak=npeak,
+					temp.dir=temp.dir,
+					chrname.rm.pattern=chrname.rm.pattern,
+					output.odir=output.odir,
+					output.npeak.file=output.npeak.file,
+					output.rpeak.file=output.rpeak.file,
+					output.rdata.file=output.rdata.file,
+					output.plot.file=output.plot.file,
+					output.result.file=output.result.file,
+					replace.flag=replace.flag,
+					clean.files.flag=clean.files.flag))
+} # end: parse.arguments()
+
+read.align <- function(align.filename) {
+# ===================================
+# Function will read a tagAlign or BAM file
+# ===================================
+	if (grepl('(\\.bam)?.*(\\.tagAlign)',align.filename)) { # if tagalign file
+		chip.data <- read.tagalign.tags(align.filename)
+		# get readlength info
+		tmpDataRows <- read.table(align.filename,nrows=500)
+		chip.data$read.length <- round(median(tmpDataRows$V3 - tmpDataRows$V2))
+	} else if (grepl('(\\.tagAlign)?.*(\\.bam)',align.filename)) { # if bam file
+		# create BAM file name
+		bam2align.filename <- sub('\\.bam','.tagAlign',align.filename)
+		# generate command to convert bam to tagalign
+		command <- vector(length=2)
+		command[1] <- sprintf("samtools view -F 0x0204 -o - %s",align.filename)
+		command[2] <- paste("awk 'BEGIN{FS=" , '"\t"' , ";OFS=", '"\t"} {if (and($2,16) > 0) {print $3,($4-1),($4-1+length($10)),"N","1000","-"} else {print $3,($4-1),($4-1+length($10)),"N","1000","+"}}', "' 1> ", bam2align.filename, sep="")
+		# command[2] <- paste("awk 'BEGIN{OFS=", '"\t"} {if (and($2,16) > 0) {print $3,($4-1),($4-1+length($10)),"N","1000","-"} else {print $3,($4-1),($4-1+length($10)),"N","1000","+"}}', "' 1> ", bam2align.filename, sep="")
+		command <- paste(command,collapse=" | ")
+		# Run command
+		status <- system(command,intern=FALSE,ignore.stderr=FALSE)
+		if ((status != 0) || !file.exists(bam2align.filename)) {
+			cat(sprintf("Error converting BAM to tagalign file: %s\n",align.filename),file=stderr())
+			q(save="no",status=1)
+		}
+		# read converted BAM file
+		chip.data <- read.tagalign.tags(bam2align.filename)
+		# get readlength info
+		tmpDataRows <- read.table(bam2align.filename,nrows=500)
+		chip.data$read.length <- round(median(tmpDataRows$V3 - tmpDataRows$V2))
+		# delete temporary tagalign file
+		file.remove(bam2align.filename)
+	} else {
+		cat(sprintf("Error:Unknown file format for file:%s\n",align.fname),file=stderr())
+		q(save="no",status=1)
+	}
+	return(chip.data)
+} # end: read.align()
+
+print.run.params <- function(params){
+# ===================================
+# Output run parameters
+# ===================================
+	cat('################\n',file=stdout())
+	cat(iparams$chip.file,
+			iparams$control.file,
+			iparams$sep.range,
+			iparams$sep.peak,
+			iparams$ex.range,
+			iparams$n.nodes,
+			iparams$fdr,
+			iparams$npeak,
+			iparams$output.odir,
+			iparams$output.npeak.file,
+			iparams$output.rpeak.file,
+			iparams$output.rdata.file,
+			iparams$output.plot.file,
+			iparams$output.result.file,
+			iparams$replace.flag,
+			labels=c('ChIP data:','Control data:', 'strandshift(min):','strandshift(step):','strandshift(max)','user-defined peak shift',
+					'exclusion(min):','exclusion(max):','num parallel nodes:','FDR threshold:','NumPeaks Threshold:','Output Directory:',
+					'narrowPeak output file name:', 'regionPeak output file name:', 'Rdata filename:',
+					'plot pdf filename:','result filename:','Overwrite files?:'),
+			fill=18,
+			file=stdout())
+	cat('\n',file=stdout())
+} # end: print.run.parameters()
+
+check.replace.flag <- function(params){
+# ===================================
+# Check if files exist
+# ===================================
+# If replace.flag is NOT set, check if output files exist and abort if necessary
+	if (! iparams$replace.flag) {
+		if (! is.na(iparams$output.npeak.file)) {
+			if (file.exists(iparams$output.npeak.file)) {
+				cat('narrowPeak file already exists. Aborting Run. Use -rf if you want to overwrite\n',file=stderr())
+				q(save="no",status=1)
+			}
+		}
+		if (! is.na(iparams$output.rpeak.file)) {
+			if (file.exists(iparams$output.rpeak.file)) {
+				cat('regionPeak file already exists. Aborting Run. Use -rf if you want to overwrite\n',file=stderr())
+				q(save="no",status=1)
+			}
+		}
+		if (! is.na(iparams$output.plot.file)) {
+			if (file.exists(iparams$output.plot.file)) {
+				cat('Plot file already exists. Aborting Run. Use -rf if you want to overwrite\n',file=stderr())
+				q(save="no",status=1)
+			}
+		}
+		if (! is.na(iparams$output.rdata.file)) {
+			if (file.exists(iparams$output.rdata.file)) {
+				cat('Rdata file already exists. Aborting Run. Use -rf if you want to overwrite\n',file=stderr())
+				q(save="no",status=1)
+			}
+		}
+	}
+}
+
+# #############################################################################
+# MAIN FUNCTION
+# #############################################################################
+
+# Check number of arguments
+minargs = 1;
+maxargs = 17;
+if (nargs < minargs | nargs > maxargs) {
+	print.usage()
+	q(save="no",status=1)
+}
+
+# Parse arguments
+# iparams$chip.file
+# iparams$isurl.chip.file
+# iparams$control.file
+# iparams$isurl.control.file
+# iparams$sep.range
+# iparams$sep.peak
+# iparams$ex.range
+# iparams$n.nodes
+# iparams$fdr
+# iparams$npeak
+# iparams$temp.dir
+# iparams$output.odir
+# iparams$output.npeak.file
+# iparams$output.rpeak.file
+# iparams$output.rdata.file
+# iparams$output.plot.file
+# iparams$output.result.file
+# iparams$replace.flag
+# iparams$clean.files.flag
+iparams <- parse.arguments(args)
+
+# Print run parameters
+print.run.params(iparams)
+
+# Check if output files exist
+check.replace.flag(iparams)
+
+# curr.chip.file and curr.control.file always point to the original ChIP and control files on disk
+# ta.chip.filename & ta.control.filename always point to the final but temporary versions of the ChIP and control files that will be passed to read.align
+
+# Download ChIP and control files if necessary to temp.dir
+if (iparams$isurl.chip.file) {
+	curr.chip.file <- file.path(iparams$temp.dir, get.file.parts(iparams$chip.file)$fullname) # file is downloaded to temp.dir. Has same name as URL suffix
+	cat('Downloading ChIP file:',iparams$chip.file,"\n",file=stdout())
+	if (system(paste('wget -N -q -P',iparams$temp.dir,iparams$chip.file)) != 0) {
+		stop('Error downloading ChIP file:',iparams$chip.file)
+	}
+} else {
+	curr.chip.file <- iparams$chip.file # file is in original directory
+}
+
+if (iparams$isurl.control.file) {
+	curr.control.file <- file.path(iparams$temp.dir, get.file.parts(iparams$control.file)$fullname) # file is downloaded to temp.dir. Has same name as URL suffix
+	cat('Downloading control file:',iparams$control.file,"\n",file=stdout())
+	if (system(paste('wget -N -q -P',iparams$temp.dir,iparams$control.file)) != 0) {
+		stop('Error downloading Control file:',iparams$control.file)
+	}
+} else {
+	curr.control.file <- iparams$control.file # file is in original directory
+}
+
+# unzip ChIP and input files if required AND copy to temp directory
+if (get.file.parts(curr.chip.file)$ext == '.gz') {
+	ta.chip.filename <- tempfile(get.file.parts(curr.chip.file)$name, tmpdir=iparams$temp.dir) # unzip file to temp.dir/[filename with .gz removed][randsuffix]
+	cat('Decompressing ChIP file\n',file=stdout())
+	if (system(paste("gunzip -c",curr.chip.file,">",ta.chip.filename)) != 0) {
+		stop('Unable to decompress file:', iparams$chip.file)
+	}
+	if (iparams$clean.files.flag) { # Remove original file if clean.files.flag is set
+		file.remove(curr.chip.file)
+	}
+} else {
+	ta.chip.filename <- tempfile(get.file.parts(curr.chip.file)$fullname, tmpdir=iparams$temp.dir)
+	if (iparams$clean.files.flag) {
+		file.rename(curr.chip.file,ta.chip.filename) # move file to temp.dir/[filename][randsuffix]
+	} else {
+		file.copy(curr.chip.file,ta.chip.filename) # copy file to temp.dir/[filename][randsuffix]
+	}
+}
+
+if (! is.na(iparams$control.file)) {
+	if (get.file.parts(curr.control.file)$ext == '.gz') {
+		ta.control.filename <- tempfile(get.file.parts(curr.control.file)$name, tmpdir=iparams$temp.dir) # unzip file to temp.dir/[filename with .gz removed][randsuffix]
+		cat('Decompressing control file\n',file=stdout())
+		if (system(paste("gunzip -c",curr.control.file,">",ta.control.filename)) != 0) {
+			stop('Unable to decompress file:', iparams$control.file)
+		}
+		if (iparams$clean.files.flag) { # Remove original file if clean.files.flag is set
+			file.remove(curr.control.file)
+		}
+	} else {
+		ta.control.filename <- tempfile(get.file.parts(curr.control.file)$fullname, tmpdir=iparams$temp.dir) # copy file to temp.dir/[filename][randsuffix]
+
+		if (iparams$clean.files.flag) {
+			file.rename(curr.control.file,ta.control.filename) # move file to temp.dir/[filename][randsuffix]
+		} else {
+			file.copy(curr.control.file,ta.control.filename) # copy file to temp.dir/[filename][randsuffix]
+		}
+	}
+}
+
+# Remove downloaded files
+if (iparams$isurl.chip.file & file.exists(curr.chip.file))  {
+	file.remove(curr.chip.file)
+}
+
+if (! is.na(iparams$control.file)) {
+	if (iparams$isurl.control.file & file.exists(curr.control.file)) {
+		file.remove(curr.control.file)
+	}
+}
+
+# Load SPP library
+library(spp)
+
+# Read ChIP tagAlign/BAM files
+cat("Reading ChIP tagAlign/BAM file",iparams$chip.file,"\n",file=stdout())
+chip.data <- read.align(ta.chip.filename)
+cat("ChIP data read length",chip.data$read.length,"\n",file=stdout())
+file.remove(ta.chip.filename) # Delete temporary file
+if (length(chip.data$tags)==0) {
+	stop('Error in ChIP file format:', iparams$chip.file)
+}
+# Remove illegal chromosome names
+if (! is.na(iparams$chrname.rm.pattern)) {
+	selectidx <- which(grepl(iparams$chrname.rm.pattern,names(chip.data$tags))==FALSE)
+	chip.data$tags <- chip.data$tags[selectidx]
+	chip.data$quality <- chip.data$quality[selectidx]
+}
+chip.data$num.tags <- sum(unlist(lapply(chip.data$tags,function(d) length(d))))
+
+# Read Control tagAlign/BAM files
+if (! is.na(iparams$control.file)) {
+	cat("Reading Control tagAlign/BAM file",iparams$control.file,"\n",file=stdout())
+	control.data <- read.align(ta.control.filename)
+	file.remove(ta.control.filename) # Delete temporary file
+	if (length(control.data$tags)==0) {
+		stop('Error in control file format:', iparams$chip.file)
+	}
+	cat("Control data read length",control.data$read.length,"\n",file=stdout())
+	# Remove illegal chromosome names
+	if (! is.na(iparams$chrname.rm.pattern)) {
+		selectidx <- which(grepl(iparams$chrname.rm.pattern,names(control.data$tags))==FALSE)
+		control.data$tags <- control.data$tags[selectidx]
+		control.data$quality <- control.data$quality[selectidx]
+	}
+	control.data$num.tags <- sum(unlist(lapply(control.data$tags,function(d) length(d))))
+}
+
+# Open multiple processes if required
+if (is.na(iparams$n.nodes)) {
+	cluster.nodes <- NULL
+} else {
+	library(snow)
+	cluster.nodes <- makeCluster(iparams$n.nodes)
+}
+
+# #################################
+# Calculate cross-correlation for various strand shifts
+# #################################
+cat("Calculating peak characteristics\n",file=stdout())
+# crosscorr
+# $cross.correlation : Cross-correlation profile as an $x/$y data.frame
+# $peak : Position ($x) and height ($y) of automatically detected cross-correlation peak.
+# $whs: Optimized window half-size for binding detection (based on the width of the cross-correlation peak)
+crosscorr <- get.binding.characteristics(chip.data,
+		srange=iparams$sep.range[c(1,3)],
+		bin=iparams$sep.range[2],
+		accept.all.tags=T,
+		cluster=cluster.nodes)
+if (!is.na(iparams$n.nodes)) {
+	stopCluster(cluster.nodes)
+}
+
+# Smooth the cross-correlation curve if required
+cc <- crosscorr$cross.correlation
+crosscorr$min.cc <- crosscorr$cross.correlation[ which.min(crosscorr$cross.correlation$y) , ] # minimum value and shift of cross-correlation
+cat("Minimum cross-correlation value", crosscorr$min.cc$y,"\n",file=stdout())
+cat("Minimum cross-correlation shift", crosscorr$min.cc$x,"\n",file=stdout())
+sbw <- 2*floor(ceiling(5/iparams$sep.range[2]) / 2) + 1 # smoothing bandwidth
+cc$y <- runmean(cc$y,sbw,alg="fast")
+
+# Compute cross-correlation peak
+bw <- ceiling(2/iparams$sep.range[2]) # crosscorr[i] is compared to crosscorr[i+/-bw] to find peaks
+peakidx <- (diff(cc$y,bw)>=0) # cc[i] > cc[i-bw]
+peakidx <- diff(peakidx,bw)
+peakidx <- which(peakidx==-1) + bw
+
+# exclude peaks from the excluded region
+if ( is.nan(iparams$ex.range[2]) ) {
+	iparams$ex.range[2] <- chip.data$read.length+10
+}
+peakidx <- peakidx[(cc$x[peakidx] < iparams$ex.range[1]) | (cc$x[peakidx] > iparams$ex.range[2])]
+cc <- cc[peakidx,]
+
+# Find max peak position and other peaks within 0.9*max_peakvalue that are further away from maxpeakposition
+maxpeakidx <- which.max(cc$y)
+maxpeakshift <- cc$x[maxpeakidx]
+maxpeakval <- cc$y[maxpeakidx]
+peakidx <-which((cc$y >= 0.9*maxpeakval) & (cc$x >= maxpeakshift))
+cc <- cc[peakidx,]
+
+# sort the peaks and get the top 3
+sortidx <- order(cc$y,decreasing=TRUE)
+sortidx <- sortidx[c(1:min(3,length(sortidx)))]
+cc.peak <- cc[sortidx,]
+
+# Override peak shift if user supplies peak shift
+if (! is.na(iparams$sep.peak)) {
+	cc.peak <- approx(crosscorr$cross.correlation$x,crosscorr$cross.correlation$y,iparams$sep.peak,rule=2)
+}
+cat("Peak cross-correlation value", paste(cc.peak$y,collapse=","),"\n",file=stdout())
+cat("Peak strand shift",paste(cc.peak$x,collapse=","),"\n",file=stdout())
+
+# Reset values in crosscorr
+crosscorr$peak$x <- cc.peak$x[1]
+crosscorr$peak$y <- cc.peak$y[1]
+
+# Compute window half size
+whs.thresh <- crosscorr$min.cc$y + (crosscorr$peak$y - crosscorr$min.cc$y)/3
+crosscorr$whs <- max(crosscorr$cross.correlation$x[crosscorr$cross.correlation$y >= whs.thresh])
+cat("Window half size",crosscorr$whs,"\n",file=stdout())
+
+# Compute phantom peak coefficient
+ph.peakidx <- which( ( crosscorr$cross.correlation$x >= ( chip.data$read.length - round(2*iparams$sep.range[2]) ) ) &
+                     ( crosscorr$cross.correlation$x <= ( chip.data$read.length + round(1.5*iparams$sep.range[2]) ) ) )
+ph.peakidx <- ph.peakidx[ which.max(crosscorr$cross.correlation$y[ph.peakidx]) ]
+crosscorr$phantom.cc <- crosscorr$cross.correlation[ph.peakidx,]
+cat("Phantom peak location",crosscorr$phantom.cc$x,"\n",file=stdout())
+cat("Phantom peak Correlation",crosscorr$phantom.cc$y,"\n",file=stdout())
+crosscorr$phantom.coeff <- crosscorr$peak$y / crosscorr$phantom.cc$y
+crosscorr$phantom.coeff <- crosscorr$peak$y / crosscorr$min.cc$y
+cat("Normalized cross-correlation coefficient (NCCC)",crosscorr$phantom.coeff,"\n",file=stdout())
+crosscorr$rel.phantom.coeff <- (crosscorr$peak$y - crosscorr$min.cc$y) / (crosscorr$phantom.cc$y - crosscorr$min.cc$y)
+cat("Relative Cross correlation Coefficient (RCCC)",crosscorr$rel.phantom.coeff,"\n",file=stdout())
+crosscorr$phantom.quality.tag <- NA
+if ( (crosscorr$rel.phantom.coeff >= 0) & (crosscorr$rel.phantom.coeff < 0.25) ) {
+	crosscorr$phantom.quality.tag <- -2
+} else if ( (crosscorr$rel.phantom.coeff >= 0.25) & (crosscorr$rel.phantom.coeff < 0.5) ) {
+	crosscorr$phantom.quality.tag <- -1
+} else if ( (crosscorr$rel.phantom.coeff >= 0.5) & (crosscorr$rel.phantom.coeff < 1) ) {
+	crosscorr$phantom.quality.tag <- 0
+} else if ( (crosscorr$rel.phantom.coeff >= 1) & (crosscorr$rel.phantom.coeff < 1.5) ) {
+	crosscorr$phantom.quality.tag <- 1
+} else if ( (crosscorr$rel.phantom.coeff >= 1.5) ) {
+	crosscorr$phantom.quality.tag <- 2
+}
+cat("Phantom Peak Quality Tag",crosscorr$phantom.quality.tag,"\n",file=stdout())
+
+# Output result to result file if required
+#Filename\tnumReads\tPeak_shift\tPeak_Correlation\tRead_length\tPhantomPeak_Correlation\tMin_Correlation_Shift\tMin_Correlation\tNormalized_CrossCorrelation_Coefficient\tRelative_CrossCorrelation_Coefficient\tQualityTag)
+if (! is.na(iparams$output.result.file)) {
+	cat(get.file.parts(iparams$chip.file)$fullname,
+			chip.data$num.tags,
+			paste(cc.peak$x,collapse=","),
+			paste(cc.peak$y,collapse=","),
+			crosscorr$phantom.cc$x,
+			crosscorr$phantom.cc$y,
+			crosscorr$min.cc$x,
+			crosscorr$min.cc$y,
+			crosscorr$phantom.coeff,
+			crosscorr$rel.phantom.coeff,
+			crosscorr$phantom.quality.tag,
+			sep="\t",
+			file=iparams$output.result.file,
+			append=TRUE)
+	cat("\n",
+			file=iparams$output.result.file,
+			append=TRUE)
+}
+
+# Save figure if required
+if (! is.na(iparams$output.plot.file)) {
+	pdf(file=iparams$output.plot.file,width=5,height=5)
+	par(mar = c(4,3.5,2,0.5), mgp = c(1.5,0.5,0), cex = 0.8);
+	plot(crosscorr$cross.correlation,
+			type='l',
+			xlab=sprintf("strand-shift (%s)",paste(cc.peak$x,collapse=",")),
+			ylab="cross-correlation")
+	abline(v=cc.peak$x,lty=2,col=2)
+	abline(v=crosscorr$phantom.cc$x,lty=2,col=4)
+	title(main=get.file.parts(iparams$chip.file)$fullname,
+	      sub=sprintf("NSC=%g,RSC=%g,Qtag=%d",crosscorr$phantom.coeff,crosscorr$rel.phantom.coeff,crosscorr$phantom.quality.tag))
+	dev.off();
+}
+
+# Save RData file if required
+if (! is.na(iparams$output.rdata.file)) {
+	save(iparams,
+			crosscorr,
+			cc.peak,
+			file=iparams$output.rdata.file);
+}
+
+# #################################
+# Call peaks
+# #################################
+
+if ( !is.na(iparams$output.npeak.file) || !is.na(iparams$output.rpeak.file) ) {
+
+	# Remove local tag anomalies
+	cat('Removing read stacks\n',file=stdout())
+	chip.data <- remove.local.tag.anomalies(chip.data$tags)
+	control.data <- remove.local.tag.anomalies(control.data$tags)
+
+	# Open multiple processes if required
+	if (is.na(iparams$n.nodes)) {
+		cluster.nodes <- NULL
+	} else {
+		cluster.nodes <- makeCluster(iparams$n.nodes)
+	}
+
+	# Find peaks
+	cat('Finding peaks\n',file=stdout())
+	if (!is.na(iparams$npeak)) {
+		iparams$fdr <- 0.96
+	}
+	narrow.peaks <- find.binding.positions(signal.data=chip.data,control.data=control.data,fdr=iparams$fdr,method=tag.lwcc,whs=crosscorr$whs,cluster=cluster.nodes)
+	if (!is.na(iparams$n.nodes)) {
+		stopCluster(cluster.nodes)
+	}
+	cat(paste("Detected",sum(unlist(lapply(narrow.peaks$npl,function(d) length(d$x)))),"peaks"),"\n",file=stdout())
+
+	# Write to narrowPeak file
+	if (!is.na(iparams$output.npeak.file)) {
+		write.narrowpeak.binding(narrow.peaks,iparams$output.npeak.file,margin=round(crosscorr$whs/2),npeaks=iparams$npeak)
+		system(paste('gzip -f ',iparams$output.npeak.file))
+	}
+
+	# Compute and write regionPeak file
+	if (!is.na(iparams$output.rpeak.file)) {
+		region.peaks <- add.broad.peak.regions(chip.data,control.data,narrow.peaks,window.size=max(50,round(crosscorr$whs/4)),z.thr=10)
+		write.narrowpeak.binding(region.peaks,iparams$output.rpeak.file,margin=round(crosscorr$whs/2),npeaks=iparams$npeak)
+		system(paste('gzip -f ',iparams$output.rpeak.file))
+	}
+
+	# Save Rdata file
+	if (! is.na(iparams$output.rdata.file)) {
+		save(iparams,
+				crosscorr,
+				cc.peak,
+				narrow.peaks,
+				region.peaks,
+				file=iparams$output.rdata.file);
+	}
+
+}
+
+