From 080af519cb46e91d0d43bed7d4df39d76856969a Mon Sep 17 00:00:00 2001 From: andrewk Date: Wed, 3 Jun 2009 03:15:45 +0000 Subject: [PATCH] Added R script and uncommented a line in recal_qual.py git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@886 348d0f76-0448-11de-a6fe-93d51630548a --- R/logistic_regression.R | 16 ++++++++++++++++ python/recal_qual.py | 2 +- 2 files changed, 17 insertions(+), 1 deletion(-) create mode 100755 R/logistic_regression.R diff --git a/R/logistic_regression.R b/R/logistic_regression.R new file mode 100755 index 000000000..604323427 --- /dev/null +++ b/R/logistic_regression.R @@ -0,0 +1,16 @@ +#!/broad/tools/apps/R-2.6.0/bin/Rscript + +args<- commandArgs(TRUE) +verbose<- TRUE + +input_fileroot <- args[1] +output_fileroot <- args[2] +dinuc <- args[3] + +currentin=paste(input_fileroot,dinuc,"csv",sep="."); +currentout=paste(output_fileroot,dinuc,"parameters",sep="."); +con <- file(currentin, "r") +con_out <- file(currentout, "w+") +data<-read.table(con,header=TRUE,sep=",") +reg<-glm(indicator~1+logitQ+I(logitQ^2)+I(logitQ^3)+I(logitQ^4)+pos+I(logitQ*pos)+I((logitQ^2)*pos)+I((logitQ^3)*pos)+I((logitQ^4)*pos)+I(pos^2)+I(logitQ*(pos^2))+I((logitQ^2)*(pos^2))+I((logitQ^3)*(pos^2))+I((logitQ^4)*(pos^2))+I(pos^3)+I(logitQ*(pos^3))+I((logitQ^2)*(pos^3))+I((logitQ^3)*(pos^3))+I((logitQ^4)*(pos^3))+I(pos^4)+I(logitQ*(pos^4))+I((logitQ^2)*(pos^4))+I((logitQ^3)*(pos^4))+I((logitQ^4)*(pos^4)),family=binomial("logit"),data=data, weights=data$count) +write(coefficients(reg),file=con_out) diff --git a/python/recal_qual.py b/python/recal_qual.py index 056719b0e..d46640fb8 100755 --- a/python/recal_qual.py +++ b/python/recal_qual.py @@ -90,7 +90,7 @@ if plot: cmd("~andrewk/covariates/plot_q_emp_stated_hst.R "+central_fileroot+".empirical_v_reported_quality.csv") cmd("~andrewk/covariates/plot_qual_diff_v_cycle_dinuc.R "+central_fileroot) -#cmd("python ~andrewk/dev/Sting/trunk/python/LogRegression.py "+central_fileroot) +cmd("python ~andrewk/dev/Sting/trunk/python/LogRegression.py "+central_fileroot) if recal_only: # Stop if we are only making recalibration tables sys.exit()