skat.R

#!/usr/bin/Rscript

library('SKAT')

cat("WARNING: this script doesn't perform a sanity check on your data and therefore it's your own responsibility to provide correct input data.\n")
args = unlist(strsplit(commandArgs(trailingOnly = TRUE)," "))
basename = strsplit(basename(args), "\\.")[[1]][1]
dat <- read.table(args[1], header=TRUE, stringsAsFactors=FALSE) 
phen = ifelse(dat$PHENOTYPE == 2, 0, 1) #Convert affected = 2 to affected = 0
geno = as.matrix(dat[,-c(1:6)])
obj <- SKAT_Null_Model(phen~1, out_type="D")
skatT <- SKAT(geno, obj, method="optimal.adj")
cat(skatT$p.value)
cat("\n\n")

warnings = names(as.vector(warnings()))
wardf <- data.frame(V1=paste("warning", 1:length(warnings), ':', sep=""), V2=warnings)
output = rbind(
	data.frame(
		V1=c("individuals", "controls", "cases", paste("pval for rho=", skatT$param$rho, ":", sep=""), "test_p:"),
		V2=c(length(phen), table(phen)[[2]], table(phen)[[1]], skatT$param$p.val.each, skatT$p.value)),
	wardf)
write.table(output, paste(basename, "_SKAT_details.txt", sep=""), col.names=F, quote=F, row.names=F, sep="\t")