[Genabel-commits] r903 - in branches/ProbABEL-pacox/v.0.1-3/examples: . test_cox_ngpreds2
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed May 2 07:30:54 CEST 2012
Author: lckarssen
Date: 2012-05-02 07:30:54 +0200 (Wed, 02 May 2012)
New Revision: 903
Added:
branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R_pacox.R
Removed:
branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R.R
Modified:
branches/ProbABEL-pacox/v.0.1-3/examples/example_cox.sh
branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R_palinear.R
Log:
example_cox.sh now exits without error when it has to. Minor changes to Cox R test scripts.
Modified: branches/ProbABEL-pacox/v.0.1-3/examples/example_cox.sh
===================================================================
--- branches/ProbABEL-pacox/v.0.1-3/examples/example_cox.sh 2012-04-29 15:23:11 UTC (rev 902)
+++ branches/ProbABEL-pacox/v.0.1-3/examples/example_cox.sh 2012-05-02 05:30:54 UTC (rev 903)
@@ -1,8 +1,42 @@
-echo "analysing Cox model"
+echo "Analysing Cox model"
-../bin/pacoxph -p coxph_data.txt -d test.mldose.2 -i test.mlinfo -m test.map -c 19 -o coxph.out.txt.dose
+outfile_base="coxph.out.txt.dose"
+outfile_dose="${outfile_base}_add.out.txt"
+../bin/pacoxph -p coxph_data.txt -d test.mldose.2 -i test.mlinfo -m test.map -c 19 -o $outfile_base
-../bin/pacoxph -p coxph_data.txt -d test.mlprob --ngpreds=2 -i test.mlinfo -m test.map -c 19 -o coxph.out.txt.prob
+outfile_base="coxph.out.txt.prob"
+outfile_prob="${outfile_base}_add.out.txt"
+../bin/pacoxph -p coxph_data.txt -d test.mlprob --ngpreds=2 -i test.mlinfo -m test.map -c 19 -o $outfile_base
+# This won't work, because differences are in second/third decimal
+#diff $outfile_dose $outfile_prob
-diff coxph.out.txt.dose_add.out.txt coxph.out.txt.prob_add.out.txt
+tmp_dose=`tempfile`
+tmp_prob=`tempfile`
+
+gawk 'NR>1 {
+ for (i=1; i<=8; i++)
+ {
+ printf "%s ", $i
+ };
+ for (i=9; i<=NF; i++)
+ {
+ printf "%6.2f ", $i
+ };
+ printf "\n"
+ }' $outfile_dose > $tmp_dose
+
+gawk 'NR>1 {
+ for (i=1; i<=8; i++)
+ {
+ printf "%s ", $i
+ };
+ for (i=9; i<=NF; i++)
+ {
+ printf "%6.2f ", $i
+ };
+ printf "\n"
+ }' $outfile_prob > $tmp_prob
+
+diff $tmp_dose $tmp_prob
+rm $tmp_dose $tmp_prob
Deleted: branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R.R
===================================================================
--- branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R.R 2012-04-29 15:23:11 UTC (rev 902)
+++ branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R.R 2012-05-02 05:30:54 UTC (rev 903)
@@ -1,92 +0,0 @@
-library(survival)
-
-###
-# load the data
-###
-
-# load phenotypic data
-pheno <- read.table("../coxph_data.txt",head=TRUE,string=FALSE)
-
-# load genetic DOSE data
-dose <- read.table("../test.mldose.2",head=FALSE,string=FALSE)
-# remove "1->" from the names of dose-IDs
-idNames <- dose[,1]
-idNames <- sub("[0-9]+->","",idNames)
-dose[,1] <- idNames
-# check consistency of names
-table(dose[,1] == pheno[,1])
-
-# load genetic PROB data
-prob <- read.table("../test.mlprob",head=FALSE,string=FALSE)
-# remove "1->" from the names of dose-IDs
-idNames <- prob[,1]
-idNames <- sub("[0-9]+->","",idNames)
-prob[,1] <- idNames
-# check consistency of names
-table(prob[,1] == pheno[,1])
-
-# check consistency DOSE <-> PROB
-doseFromProb <- matrix(NA,ncol=dim(dose)[2],nrow=dim(dose)[1])
-for (i in 3:dim(dose)[2]) {
- indexHom <- 3 + ( i - 3 ) * 2
- indexHet <- indexHom + 1
- doseFromProb[,i] <- prob[,indexHom]*2 + prob[,indexHet]
- print( table( ( dose[,i] - doseFromProb[,i] ) > 1e-8 ) )
-}
-
-###
-# run analysis
-###
-
-# run ProbABEL
-system("cd ..; sh example_cox.sh; cd test_cox_ngpreds2")
-resPaAddDose <- read.table("../coxph.out.txt.dose_add.out.txt", head=TRUE)[,c("beta_SNP_add","sebeta_SNP_add")]
-resPaAddProb <- read.table("../coxph.out.txt.prob_add.out.txt", head=TRUE)[,c("beta_SNP_addA1","sebeta_SNP_addA1")]
-resPaDom <- read.table("../coxph.out.txt.prob_domin.out.txt", head=TRUE)[,c("beta_SNP_domA1","sebeta_SNP_domA1")]
-resPaRec <- read.table("../coxph.out.txt.prob_recess.out.txt", head=TRUE)[,c("beta_SNP_recA1","sebeta_SNP_recA1")]
-resPaOdom <- read.table("../coxph.out.txt.prob_over_domin.out.txt", head=TRUE)[,c("beta_SNP_odom","sebeta_SNP_odom")]
-
-attach(pheno)
-
-# run analysis on dose
-emptyRes <- resPaAddDose
-emptyRes[] <- NA
-resRAddDose <- resRAddProb <- resRDom <- resRRec <- resROdom <- emptyRes
-for (i in 3:dim(dose)[2]) {
- smr <- summary( coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + dose[,i] ) )
- resRAddDose[i-2,] <- smr$coef[4,c(1,3)]
- smr <- summary( coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + doseFromProb[,i] ) )
- resRAddProb[i-2,] <- smr$coef[4,c(1,3)]
-}
-resPaAddDose/resRAddDose
-resPaAddProb/resRAddProb
-
-# dominant model
-for (i in 3:dim(dose)[2]) {
- indexHom <- 3 + ( i - 3 ) * 2
- indexHet <- indexHom + 1
- regProb <- prob[,indexHom] + prob[,indexHet]
- smr <- summary( coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + regProb ) )
- resRDom[i-2,] <- smr$coef[4,c(1,3)]
-}
-resPaDom/resRDom
-
-# recessive model
-for (i in 3:dim(dose)[2]) {
- indexHom <- 3 + ( i - 3 ) * 2
- indexHet <- indexHom + 1
- regProb <- prob[,indexHom]
- smr <- summary( coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + regProb ) )
- resRRec[i-2,] <- smr$coef[4,c(1,3)]
-}
-resPaRec/resRRec
-
-# over-dominant model
-for (i in 3:dim(dose)[2]) {
- indexHom <- 3 + ( i - 3 ) * 2
- indexHet <- indexHom + 1
- regProb <- prob[,indexHet]
- smr <- summary( coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + regProb ) )
- resROdom[i-2,] <- smr$coef[4,c(1,3)]
-}
-resPaOdom/resROdom
Copied: branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R_pacox.R (from rev 902, branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R.R)
===================================================================
--- branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R_pacox.R (rev 0)
+++ branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R_pacox.R 2012-05-02 05:30:54 UTC (rev 903)
@@ -0,0 +1,92 @@
+library(survival)
+
+###
+# load the data
+###
+
+# load phenotypic data
+pheno <- read.table("../coxph_data.txt",head=TRUE,string=FALSE)
+
+# load genetic DOSE data
+dose <- read.table("../test.mldose.2",head=FALSE,string=FALSE)
+# remove "1->" from the names of dose-IDs
+idNames <- dose[,1]
+idNames <- sub("[0-9]+->","",idNames)
+dose[,1] <- idNames
+# check consistency of names
+table(dose[,1] == pheno[,1])
+
+# load genetic PROB data
+prob <- read.table("../test.mlprob",head=FALSE,string=FALSE)
+# remove "1->" from the names of prob-IDs
+idNames <- prob[,1]
+idNames <- sub("[0-9]+->","",idNames)
+prob[,1] <- idNames
+# check consistency of names
+table(prob[,1] == pheno[,1])
+
+# check consistency DOSE <-> PROB
+doseFromProb <- matrix(NA, ncol=dim(dose)[2], nrow=dim(dose)[1])
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ doseFromProb[,i] <- prob[,indexHom]*2 + prob[,indexHet]
+ print( table( ( dose[,i] - doseFromProb[,i] ) > 1e-8 ) )
+}
+
+###
+# run analysis
+###
+
+# run ProbABEL
+system("cd ..; sh example_cox.sh; cd test_cox_ngpreds2")
+resPaAddDose <- read.table("../coxph.out.txt.dose_add.out.txt", head=TRUE)[,c("beta_SNP_add","sebeta_SNP_add")]
+resPaAddProb <- read.table("../coxph.out.txt.prob_add.out.txt", head=TRUE)[,c("beta_SNP_addA1","sebeta_SNP_addA1")]
+resPaDom <- read.table("../coxph.out.txt.prob_domin.out.txt", head=TRUE)[,c("beta_SNP_domA1","sebeta_SNP_domA1")]
+resPaRec <- read.table("../coxph.out.txt.prob_recess.out.txt", head=TRUE)[,c("beta_SNP_recA1","sebeta_SNP_recA1")]
+resPaOdom <- read.table("../coxph.out.txt.prob_over_domin.out.txt", head=TRUE)[,c("beta_SNP_odom","sebeta_SNP_odom")]
+
+attach(pheno)
+
+# run analysis on dose
+emptyRes <- resPaAddDose
+emptyRes[] <- NA
+resRAddDose <- resRAddProb <- resRDom <- resRRec <- resROdom <- emptyRes
+for (i in 3:dim(dose)[2]) {
+ smr <- summary( coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + dose[,i] ) )
+ resRAddDose[i-2,] <- smr$coef[4,c(1,3)]
+ smr <- summary( coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + doseFromProb[,i] ) )
+ resRAddProb[i-2,] <- smr$coef[4,c(1,3)]
+}
+resPaAddDose/resRAddDose
+resPaAddProb/resRAddProb
+
+# dominant model
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ regProb <- prob[,indexHom] + prob[,indexHet]
+ smr <- summary( coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + regProb ) )
+ resRDom[i-2,] <- smr$coef[4,c(1,3)]
+}
+resPaDom/resRDom
+
+# recessive model
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ regProb <- prob[,indexHom]
+ smr <- summary( coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + regProb ) )
+ resRRec[i-2,] <- smr$coef[4,c(1,3)]
+}
+resPaRec/resRRec
+
+# over-dominant model
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ regProb <- prob[,indexHet]
+ smr <- summary( coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + regProb ) )
+ resROdom[i-2,] <- smr$coef[4,c(1,3)]
+}
+resPaOdom/resROdom
Modified: branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R_palinear.R
===================================================================
--- branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R_palinear.R 2012-04-29 15:23:11 UTC (rev 902)
+++ branches/ProbABEL-pacox/v.0.1-3/examples/test_cox_ngpreds2/run_models_in_R_palinear.R 2012-05-02 05:30:54 UTC (rev 903)
@@ -16,7 +16,7 @@
# load genetic PROB data
prob <- read.table("../test.mlprob",head=FALSE,string=FALSE)
-# remove "1->" from the names of dose-IDs
+# remove "1->" from the names of prob-IDs
idNames <- prob[,1]
idNames <- sub("[0-9]+->","",idNames)
prob[,1] <- idNames
@@ -24,7 +24,7 @@
table(prob[,1] == pheno[,1])
# check consistency DOSE <-> PROB
-doseFromProb <- matrix(NA,ncol=dim(dose)[2],nrow=dim(dose)[1])
+doseFromProb <- matrix(NA, ncol=dim(dose)[2], nrow=dim(dose)[1])
for (i in 3:dim(dose)[2]) {
indexHom <- 3 + ( i - 3 ) * 2
indexHet <- indexHom + 1
More information about the Genabel-commits
mailing list