diff xcorrhelper.r @ 0:d42f4d78c85e draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/xcorrhelper.r	Wed Dec 17 10:40:23 2014 -0500
@@ -0,0 +1,111 @@
+#######################################
+#### Processes and saves xcorr output
+#### R script		
+######################################
+
+
+
+# ============================================== 
+# Parsing Arguments (source: phantom SPP script)
+# ============================================== 
+args = commandArgs(trailingOnly = TRUE) # Read Arguments from command line
+#nargs = length(args) # number of arguments
+
+#Set arguments to default values
+infile = NA # input file
+out = NA # output directory
+
+#Parsing arguments and storing values
+for (each.arg in args) {
+	#input bed file
+	if (grepl('^-infile=',each.arg)) {
+		arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+		if (! is.na(arg.split[2]) ) {
+				infile <- arg.split[2]
+		} else {
+			stop('No input file')
+		} 
+	}
+	if (grepl('^-out=',each.arg)) {			
+		arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
+		if (! is.na(arg.split[2]) ) {
+			out <- arg.split[2] 
+		} else {
+			stop('No output directory')
+		}
+	}
+}
+
+#=======================> DONE! 
+
+
+# =================== 
+# Do the processing
+# ===================
+
+#Read in variables
+if (file.exists(infile)) {
+
+message("@", infile)
+xcorrs = read.table(infile, header=FALSE)
+shifts = as.numeric(xcorrs$V2) 
+
+
+#how many chromosomes
+num = length(xcorrs$V1)
+
+if (length(xcorrs$V2) > 0) {
+if (length(xcorrs$V2) > 1)  {
+#shift stats
+avg = mean(shifts)
+std = sd(shifts)
+mod = as.numeric(names(sort(-table(shifts)))[1])
+
+#find out fragment length
+fraglength = 100 #default
+default = 1 #was default shift assumed
+if (length(shifts) != 0) {
+	fraglength = ceiling(mean(shifts))
+	default = 0
+} 
+} else {
+avg = xcorrs$V2
+std = xcorrs$V2
+mod = xcorrs$V2
+fraglength = xcorrs$V2
+default = 0
+}
+} else {
+fraglength = 100 #default
+}
+} else {
+fraglength = 100
+default = 1
+}
+#=======================> DONE! 
+
+
+
+
+
+# ===================== 
+# Write result to File
+# =====================
+if (default == 0) {
+writethis = paste0("Number of chromosomes used:", num, "\nAverage strand shift:", avg, "\nStandard Deviation:", std, "\nMode:", mod, "\nFragment Length:", fraglength)
+write(writethis, file = paste0(out, "/xcorrsummary.txt"))
+
+if (length(xcorrs$V2) > 1)  {
+	i = pdf(file=paste0(out, "/xcorrsummary.pdf"))
+	i = plot(density(shifts), main="Strand Shifts Distribution over Chromosomes")
+	i = dev.off()
+}
+	message("Fragment Length:", fraglength, "\n")
+} else if (default == 1) {
+
+writethis = paste0("Number of chromosomes used:", "0", "\nAverage strand shift:", "NA", "\nStandard Deviation:", "NA", "\nMode:", "NA", "\nFragment Length:", fraglength)
+write(writethis, file = paste0(out, "/xcorrsummary.txt"))
+
+	message("Fragment length estimation failed, default fragment length assumed:", fraglength, "\n")
+}
+#=======================> DONE!