zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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
This commit is contained in:
andrewk 2009-06-03 03:15:45 +00:00
parent b2eb724456
commit 080af519cb
2 changed files with 17 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
R/logistic_regression.R 100755
View File

@ -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)

2
python/recal_qual.py
View File

@ -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()