Clean R code

corescripts/generate_rsb.R

@@ -1,3 +1,4 @@
+#!/bin/env Rscript
 #
 #  Murray Cadzow and James Boocock
 #  November 2013
@@ -9,22 +10,26 @@ require(getopt)
 
 require(rehh)
 
-args=commandArgs(TRUE)
-spec = matrix(c(
-	'pop1', 	'p', 1, 'character',
-	'pop2',   'P', 1, 'character',
-	'chr',    'c', 1, 'character',
-	'pop1file','-i',1 ,'character',
-	'pop2file,','-I',1,'character'
-),byrow=T,ncol=4)
-opt = getopt(spec)	
-
-if(!is.null(opt$help)){
-	cat(getopt(spec,usage=TRUE));
-	q(status=1);
+args <- commandArgs(TRUE)
+spec <- matrix(c(
+	'pop1', 	   'p', 1, 'character',
+	'pop2',      'P', 1, 'character',
+	'chr',       'c', 1, 'character',
+	'pop1file', '-i', 1, 'character',
+	'pop2file', '-I', 1, 'character'
+               ),byrow=T,ncol=4)
+
+
+opt <- getopt(spec)
+
+if (!is.null(opt$help)) {
+  cat(getopt(spec, usage = TRUE))
+  quit(status = 1)
 }
-pop1_data = read.table(opt$pop1file,header=T,row.names=1)
-pop2_data = read.table(opt$pop2file,header=T,row.names=1)
 
-rsb_out = ies2rsb(pop1_data,pop2_data,popname1=opt$pop1,popname2=opt$pop2)
-write.table(rsb_out$res.rsb,paste(opt$chr,opt$pop1,opt$pop2,'.rsb',sep=""))
+pop1_data <- read.table(opt$pop1file, header = T, row.names = 1)
+pop2_data <- read.table(opt$pop2file, header = T, row.names = 1)
+
+rsb_out <- ies2rsb(pop1_data, pop2_data, popname1 = opt$pop1, popname2 = opt$pop2)
+
+write.table(rsb_out$res.rsb, paste(opt$chr, opt$pop1, opt$pop2, ".rsb", sep = "")) 

corescripts/haps_indel_and_maf_filter.R

@@ -1,25 +1,32 @@
+#!/bin/env Rscript
 #
 # Murray Cadzow
 # July 2013
 # University of Otago
 
-args<-commandArgs(TRUE)
-#read in haps file from shapeit
-hapsPop=read.table(file=args[1])
-maf=as.numeric(args[2])
-output=as.character(args[3])
+args <- commandArgs(TRUE)
 
-hapsPop=read.table(file=args[1])
-hapsPop=hapsPop[nchar(as.character(hapsPop[,4]))==1 &
-                     as.character(hapsPop[,4]) != "-" &
-                     nchar(as.character(hapsPop[,5]))==1 &
-                     as.character(hapsPop[,5]) != "-", ] #remove indels
-hapsPop[,1]= hapsPop[,2]
+## read in haps file from shapeit
 
+hapsPop <- read.table(file = args[1])
+maf     <- as.numeric(args[2])
+output  <- as.character(args[3])
 
-af = apply(hapsPop[,6:length(hapsPop[1,])], 1, FUN=function(x){sum(x)/length(x)}  )
-hapsPop=hapsPop[(af >= maf) & (af <= (1-maf)),]
-print(paste("af >",maf))
-print(table((af > maf ) & (af < (1-maf))))
+hapsPop <- read.table(file=args[1])
+hapsPop <- hapsPop[nchar(as.character(hapsPop[,4]))==1 &
+                   as.character(hapsPop[,4]) != "-"    &
+                   nchar(as.character(hapsPop[,5]))==1 &
+                   as.character(hapsPop[,5]) != "-", ]   # remove indels
 
-write.table(hapsPop, file=output, quote=FALSE, row.names=FALSE, col.names=FALSE)
+hapsPop[,1] <- hapsPop[,2]
+
+## This could be replaced with apply(..., FUN=mean), or rowMeans(x)
+af <- apply(hapsPop[, 6:length(hapsPop[1, ])], 1, FUN = function(x) {
+  sum(x)/length(x)
+})
+
+hapsPop <- hapsPop[(af >= maf) & (af <= (1 - maf)), ]
+print(paste("af >", maf))
+print(table((af > maf) & (af < (1 - maf))))
+
+write.table(hapsPop, file = output, quote = FALSE, row.names = FALSE, col.names = FALSE) 

corescripts/make_plots.R

@@ -1,40 +1,42 @@
+#!/bin/env Rscript
+
 require(getopt)
 
 spec = matrix(c(
-    'help',     'h',    0,  "logical",
-    'plot_type','p',    1,  "character",
-    'data_file', 'd',   1,  "character",
+    'help',             'h', 0,  "logical",
+    'plot_type',        'p', 1,  "character",
+    'data_file',        'd', 1,  "character",
     'plot_output_name', 'n', 1, 'character'
- ), byrow=T, ncol=4)
+             ), byrow=T, ncol=4)
 
-opt = getopt(spec)
+opt <- getopt(spec)
 
-if (!is.null(opt$help)){
-    cat(getopt(spec,usage=TRUE));
-    q(status=1);
+if (!is.null(opt$help)) {
+  cat(getopt(spec, usage = TRUE))
+  quit(status = 1)
 }
 
-if (opt$plot_type == "taj"){
-    tajimaD=read.table(file=opt$data_file, header=TRUE)
-    png(opt$plot_output_name)
-    plot(tajimaD[,4] ~ tajimaD[,2],pch='.',cex=2,xlab="Chromosome position (bp)", ylab="D statistic")
-    dev.off()
-}else if( opt$plot_type == "ihs"){
-    ihs = read.table(file=opt$data_file)
-    png(opt$plot_output_name)
-    plot(ihs[,4] ~ ihs[,2],,pch='.',cex=2,ylab=expression("-" * log[10] * "[" ~ "1-2|" * Phi[scriptstyle(italic(iHS))] * "-0.5|" ~ "]"),
-         xlab="Chromosome Position BP")
-    dev.off()
-}else if(opt$plot_type == "fay"){
-    fay = read.table(file=opt$data_file,comment.char="#") 
-    png(opt$plot_output_name)
-    plot(CEUFay[,15] ~ CEUFay[,1],xlab='Chromosome position (bp)',ylab="H Statistic"))
-    dev.off()
-}else if(opt$plot_type == "rsb"){
-
-}else if(opt$plot_type == "fst_weir"){
-
-}else if(opt$plot_type == "fst_hap"){
-
+if (opt$plot_type == "taj") {
+  tajimaD <- read.table(file = opt$data_file, header = TRUE)
+  png(opt$plot_output_name)
+  plot(tajimaD[, 4] ~ tajimaD[, 2], pch = ".", cex = 2, xlab = "Chromosome position (bp)", ylab = "D statistic")
+  dev.off()
+} else if (opt$plot_type == "ihs") {
+  ihs <- read.table(file = opt$data_file)
+  png(opt$plot_output_name)
+  plot(ihs[, 4] ~ ihs[, 2], , pch = ".", cex = 2, ylab = expression("-" * log[10] * "[" ~ "1-2|" * Phi[scriptstyle(italic(iHS))] * 
+                                                                      "-0.5|" ~ "]"), xlab = "Chromosome Position BP")
+  dev.off()
+} else if (opt$plot_type == "fay") {
+  fay <- read.table(file = opt$data_file, comment.char = "#")
+  png(opt$plot_output_name)
+  plot(CEUFay[, 15] ~ CEUFay[, 1], xlab = "Chromosome position (bp)", ylab = "H Statistic")
+  dev.off()
+} else if (opt$plot_type == "rsb") {
+  
+} else if (opt$plot_type == "fst_weir") {
+  
+} else if (opt$plot_type == "fst_hap") {
+  
 }