From 1e43ca81364c99eef919ad0fbf35657d46e903dd Mon Sep 17 00:00:00 2001
From: Tanushree Tunstall <tanu@tunstall.in>
Date: Tue, 23 Jun 2020 11:57:51 +0100
Subject: [PATCH] tody scracth script for various  OR calcs

---
 scripts/af_or_calcs_scratch.R | 288 ++++++++++++----------------------
 1 file changed, 104 insertions(+), 184 deletions(-)

diff --git a/scripts/af_or_calcs_scratch.R b/scripts/af_or_calcs_scratch.R
index beda144..a4b02e2 100644
--- a/scripts/af_or_calcs_scratch.R
+++ b/scripts/af_or_calcs_scratch.R
@@ -154,14 +154,19 @@ cat(paste0('Total no. of distinct comp snps to perform OR calcs: ', length(gene_
 #3) modified chisq.test
 #4) logistic
 #5) adjusted logistic?
-#6) kinship (separate script)
 
 #======================================
+#########################
+# custom chisq function: 
+# To calculate OR
+#########################
+#i = "pnca_p.trp68gly"
+#mut = grepl(i,raw_data$all_muts_gene)
+#dst = raw_data[[drug]] 
 
-################# modified chisq OR
-# Define OR function
 #x = as.numeric(mut)
 #y = dst
+
 custom_chisq_or = function(x,y){
   tab = as.matrix(table(x,y))
   a = tab[2,2]
@@ -176,9 +181,10 @@ custom_chisq_or = function(x,y){
   
 }
 
-#========================
+#======================================
+#==============
 # TEST WITH ONE
-#========================
+#==============
 i = "pnca_p.trp68gly"
 i = "pnca_p.gln10pro"
 i = "pnca_p.leu159arg"
@@ -207,84 +213,93 @@ table(sid)
 # 3X2 table
 table(mut, dst, sid)
 
-#============================
-# compare OR
+#===================================================
+# compare ORs from different calcs
+
+#1) chisq
 chisq.test(table(mut,dst))
-chisq.test(table(mut,dst)) $ statistic
+chisq_estimate = chisq.test(table(mut,dst))$statistic
+est_chisq =  chisq_estimate[[1]]; print(paste0('chi sq estimate:', est_chisq))# numeric part only 
+pval_chisq = chisq.test(table(mut,dst))$p.value; print(paste0('pvalue:', pval_chisq))
 
-f = chisq.test(table(mut,dst)) $ statistic
-chisq.test(dst, mut) $ statistic
 
+#2) fisher
 fisher.test(table(mut, dst))
 fisher.test(table(mut, dst))$p.value
-fisher.test(table(mut, dst))$estimate
-logistic_chisq_or(mut,dst)
+est_fisher = fisher.test(table(mut, dst))$estimate
+or_fisher = est_fisher[[1]]; print(paste0('OR fisher:', or_fisher))# numeric part only
+pval_fisher = fisher.test(table(mut, dst))$p.value; print(paste0('pval fisher:', pval_fisher))
 
-# logistic or
+
+#3) custom chisq
+or_mychisq = custom_chisq_or(mut,dst)
+
+#4) logistic
 summary(model<-glm(dst ~ mut, family = binomial))
-or_logistic = exp(summary(model)$coefficients[2,1]); print(or_logistic)
-pval_logistic_maxit = summary(model)$coefficients[2,4]; print(pval_logistic_maxit)
+
+or_logistic = exp(summary(model)$coefficients[2,1]); print(paste0('OR logistic:', or_logistic))
+pval_logistic = summary(model)$coefficients[2,4]; print(paste0('pval logistic:', pval_logistic))
 
 # extract SE of the logistic model for a given snp
-logistic_se = summary(model)$coefficients[2,2]
-print(paste0('SE:', logistic_se))
+logistic_se = summary(model)$coefficients[2,2]; print(paste0('SE:', logistic_se))
     
 # extract Z of the logistic model for a given snp
-logistic_zval = summary(model)$coefficients[2,3]
-print(paste0('Z-value:', logistic_zval))
-  
-  
+logistic_zval = summary(model)$coefficients[2,3]; print(paste0('Z-value:', logistic_zval))
+
 # extract confint interval of snp (2 steps, since the output is a named number)
-ci_mod = exp(confint(model))[2,]
-print(paste0('CI:', ci_mod))
+ci_mod = exp(confint(model))[2,]; print(paste0('CI:', ci_mod))
 #logistic_ci = paste(ci_mod[["2.5 %"]], ",", ci_mod[["97.5 %"]])
   
-logistic_ci_lower = ci_mod[["2.5 %"]]
-logistic_ci_upper = ci_mod[["97.5 %"]]
+logistic_ci_lower = ci_mod[["2.5 %"]]; print(paste0('CI_lower:', logistic_ci_lower))
+logistic_ci_upper = ci_mod[["97.5 %"]]; print(paste0('CI_upper:', logistic_ci_upper))
   
-print(paste0('CI_lower:', logistic_ci_lower))
-print(paste0('CI_upper:', logistic_ci_upper))
-
-
 # adjusted logistic or: doesn't seem to make a difference
 summary(model2<-glm(dst ~ mut + sid, family = binomial))
-or_logistic2 = exp(summary(model2)$coefficients[2,1]); print(or_logistic2)
-pval_logistic2 = summary(model2)$coefficients[2,4]; print(pval_logistic2)
+or_logistic2 = exp(summary(model2)$coefficients[2,1]); print(paste0('Adjusted OR logistic:', or_logistic2))
+pval_logistic2 = summary(model2)$coefficients[2,4]; print(paste0('Adjusted pval logistic:',pval_logistic2))
 
-#============ looping with sapply
-#####################
-# iterate: subset
-#####################
-
-snps_test = c("pnca_p.trp68gly", "pnca_p.leu4ser", "pnca_p.leu159arg","pnca_p.his57arg" )
-
-snps = snps_test[1:4]
+# print all output 
+writeLines(c(paste0("mutation:", i)
+             , paste0("=========================")
+             , paste0("OR logistic:", or_logistic,"--->", "pval logistic:", pval_logistic )
+             , paste0("OR adjusted logistic:", or_logistic2,"--->", "pval adjusted logistic:", pval_logistic2)
+             , paste0("OR fisher:", or_fisher, "--->","pval fisher:", pval_fisher )
+             , paste0("OR custom chisq:", or_mychisq, "--->","pval fisher:", pval_fisher )
+             , paste0("Chisq estimate:", est_chisq, "--->","pval chisq:", pval_chisq)))
+#%%========================================================
+# looping with sapply: a subset of mutations
+#%%========================================================
 
+snps = c("pnca_p.trp68gly", "pnca_p.leu4ser", "pnca_p.leu159arg","pnca_p.his57arg" )
+#snps = snps_test[1:2]
 snps
 
-
+# custom chisq
 ors = sapply(snps,function(m){
   mut = grepl(m,raw_data$all_muts_gene)
-  logistic_chisq_or(mut,dst)
+  custom_chisq_or(mut,dst)
 })
+head(ors)
 
-ors
-
+# pvalue fisher, to be used with custom chisq
 pvals = sapply(snps,function(m){
   mut = grepl(m,raw_data$all_muts_gene)
   fisher.test(mut,dst)$p.value
 })
+head(pvals)
 
-pvals
-
+# allele frequency
 afs = sapply(snps,function(m){
   mut = grepl(m,raw_data$all_muts_gene)
   mean(mut)
 })
+head(afs)
 
-afs
+# logistic reg parameters: individual sapply
 
-## logistic or
+#--------------
+## logistci or
+#--------------
  ors_logistic = sapply(snps,function(m){
   mut = grepl(m,raw_data$all_muts_gene)
   model<-glm(dst ~ mut, family = binomial)
@@ -294,7 +309,9 @@ afs
 ors_logistic
 head(ors_logistic); head(names(ors_logistic))
 
+#-------------------
 ## logistic p-value
+#--------------
 pvals_logistic = sapply(snps,function(m){
   mut = grepl(m,raw_data$all_muts_gene)
   model<-glm(dst ~ mut , family = binomial)
@@ -303,7 +320,9 @@ pvals_logistic = sapply(snps,function(m){
 
 head(pvals_logistic); head(names(pvals_logistic))
 
+#--------------
 ## logistic se
+#--------------
 se_logistic = sapply(snps,function(m){
   mut = grepl(m,raw_data$all_muts_gene)
   model<-glm(dst ~ mut , family = binomial)
@@ -312,8 +331,9 @@ se_logistic = sapply(snps,function(m){
 
 head(se_logistic); head(names(se_logistic))
 
-
+#--------------
 ## logistic z-value
+#--------------
 zval_logistic = sapply(gene_snps_unique,function(m){
   mut = grepl(m,raw_data$all_muts_gene)
   model<-glm(dst ~ mut , family = binomial)
@@ -321,8 +341,9 @@ zval_logistic = sapply(gene_snps_unique,function(m){
 })
 
 head(zval_logistic); head(names(zval_logistic))
-
+#--------------
 ## logistic ci - lower bound
+#--------------
 ci_lb_logistic = sapply(snps,function(m){
   mut = grepl(m,raw_data$all_muts_gene)
   model<-glm(dst ~ mut , family = binomial)
@@ -331,8 +352,9 @@ ci_lb_logistic = sapply(snps,function(m){
 })
 
 head(ci_lb_logistic); head(names(ci_lb_logistic))
-
+#--------------
 ## logistic ci - upper bound
+#--------------
 ci_ub_logistic = sapply(snps,function(m){
   mut = grepl(m,raw_data$all_muts_gene)
   model<-glm(dst ~ mut , family = binomial)
@@ -343,7 +365,9 @@ ci_ub_logistic = sapply(snps,function(m){
 
 head(ci_ub_logistic); head(names(ci_ub_logistic))
 
-# logistic adj # Doesn't seem to make a difference
+#--------------
+# adjusted logistic with sample id: Doesn't seem to make a difference
+#--------------
 logistic_ors2 = sapply(snps,function(m){
   mut = grepl(m,raw_data$all_muts_gene)
   c = raw_data$id[mut]
@@ -352,38 +376,34 @@ logistic_ors2 = sapply(snps,function(m){
   or_logistic2 = exp(summary(model2)$coefficients[2,1])
   #pval_logistic2 = summary(model2)$coefficients[2,4]
 })
-
-logistic_ors2
-
-or_logistic2; pval_logistic2
-
 head(logistic_ors)
-#===========================================================
-#%%
-# sapply with multiple values
 
-#https://gist.github.com/primaryobjects/33adabc337edd67b4a8d
-
-snps_test = c("pnca_p.trp68gly", "pnca_p.leu4ser", "pnca_p.leu159arg","pnca_p.his57arg" )
-
-snps = snps_test[1:4]
+#=!=!=!=!=!=!=!=!=!=!=!=!=!=!
+# COMMENT: individual sapply seem wasteful
+#=!=!=!=!=!=!=!=!=!=!=!=!=!=!
 
+#%%========================================================
+# sapply with multiple values being returned as df
+#Link: https://gist.github.com/primaryobjects/33adabc337edd67b4a8d
+#%%========================================================
+snps = c("pnca_p.trp68gly", "pnca_p.leu4ser", "pnca_p.leu159arg","pnca_p.his57arg" )
+#snps = snps_test[1:4]
 snps
 
+# yayy works!
 # DV: pyrazinamide 0 or 1
 dst = raw_data[[drug]]
 
-# yayy works!
-testdf = data.frame()
+# initialise an empty df
+or_df = data.frame()
 
 x = sapply(snps,function(m){
 
-  df = data.frame()
+  #df = data.frame()
   
   mut = grepl(m,raw_data$all_muts_gene)
   model<-glm(dst ~ mut, family = binomial)
   
-  
   # allele frequency
   afs = mean(mut)
   
@@ -419,7 +439,7 @@ x = sapply(snps,function(m){
   
   pval_chisq = chisq.test(dst, mut)$p.value                                                        
                                                          
-  #build a row to append to df
+  # build a row to append to df
   row = data.frame(mutation = m
                    , af = afs
                    , beta_logistic = beta_logistic
@@ -439,15 +459,15 @@ x = sapply(snps,function(m){
                    )
   #print(row)
   
-  testdf <<- rbind(testdf, row)
+  or_df <<- rbind(or_df, row)
   
 })
 
-write.csv(testdf, 'test_ors.csv')
+write.csv(or_df, 'test_ors.csv')
+
 #=================================
-####################
-# iterate: subset
-#####################
+# testing logistic or with maxit, etc.
+
 print(paste0('subset to iterate over;', snps))
 
 # start loop
@@ -499,7 +519,12 @@ for (i in snps){
   or_fisher = or_fisher[[1]]; or_fisher
   
   pval_fisher = fisher.test(table(dst, mut))$p.value ; pval_fisher
-  
+
+  #=====================
+  # custom chi square
+  #=====================
+  or_mychisq = custom_chisq_or(mut,dst)
+
   #=====================
   # chi square
   #=====================
@@ -513,9 +538,10 @@ for (i in snps){
   # all output 
   writeLines(c(paste0("mutation:", i)
                , paste0("=========================")
-               , paste0("or_logistic:", or_logistic,"--->", "P-val_logistic_maxit:", pval_logistic_maxit )
-               , paste0("OR_fisher:", or_fisher, "--->","P-val_fisher:", pval_fisher )
-               , paste0("Chi_sq_estimate:", est_chisq, "--->","P-val_chisq:", pval_chisq)))
+               , paste0("OR logistic:", or_logistic,"--->", "pval logistic:", pval_logistic )
+               , paste0("OR fisher:", or_fisher, "--->","pval fisher:", pval_fisher )
+               , paste0("OR custom chisq:", or_mychisq, "--->","pval fisher:", pval_fisher )
+               , paste0("Chisq estimate:", est_chisq, "--->","pval chisq:", pval_chisq)))
 }
   #=====================
   # fishers test
@@ -529,110 +555,4 @@ for (i in snps){
   
  # https://stats.stackexchange.com/questions/259635/what-is-the-difference-using-a-fishers-exact-test-vs-a-logistic-regression-for
   exact2x2(table(dst, mut),tsmethod="central")
-
-
 #=====================================================================
-# iterate over a df and then add these values  
-#
-my_data = as.data.frame(gene_snps_unique)
-colnames(my_data) = "mutation"
-print(colnames(my_data))
-
-perfectSeparation <- function(w) {
-  if(grepl("fitted probabilities numerically 0 or 1 occurred",
-           as.character(w))) {ww <<- ww+1}
-}
-
-for(i in my_data$mutation) {
-  print(paste0('snp to iterate over:', i))
-}
-
-for(i in my_data$mutation) {
-  print(paste0('snp to iterate over:', i))
-
-  #####
-  # Run logistic regression
-  #####
-  
-  #*************
-  # start logistic regression model building
-  # set the IV and DV for the logistic regression model and model
-  #*************
-  # IV: corresponds to each unique snp (extracted using grep)
-  mut = as.numeric(grepl(i,raw_data$dr_muts_pza))
-  
-  # DV: pyrazinamide 0 or 1
-  dst = as.numeric(raw_data$pyrazinamide)
-  
-  tab = table(mut, dst)
-  print(tab)
-  
-  # glm model: with and without maxit
-  model = tryCatch( glm(dst ~ mut
-                        , family = binomial
-                        #, control = glm.control(maxit = 1) # only used when required for one step estimator
-  ), warning = perfectSeparation)
-  
-  model = glm(dst ~ mut, family = binomial)
-  
-  print(summary(model))
-  
-  #**********
-  # extract relevant model output
-  #**********
-  # extract log OR i.e the Beta estimate of the logistic model for a given snp
-  my_logor = summary(model)$coefficients[2,1]
-  print(paste0('Beta:', my_logor))
-  
-  # Dervive OR i.e exp(my_or) from the logistic model for a given snp
-  #my_or = round(exp(summary(model)$coefficients[2,1]), roundto)
-  my_or = exp(summary(model)$coefficients[2,1])
-  print(paste0('OR:', my_or))
-  
-  # extract SE of the logistic model for a given snp
-  my_se = summary(model)$coefficients[2,2]
-  print(paste0('SE:', my_se))
-  
-  # extract Z of the logistic model for a given snp
-  my_zval = summary(model)$coefficients[2,3]
-  print(paste0('Z-value:', my_zval))
-
-  # extract P-value of the logistic model for a given snp
-  my_pval = summary(model)$coefficients[2,4]
-  print(paste0('P-value:', my_pval))
-  
-  # extract confint interval of snp (2 steps, since the output is a named number)
-  ci_mod = exp(confint(model))[2,]
-  #my_ci = paste(ci_mod[["2.5 %"]], ",", ci_mod[["97.5 %"]])
-  
-  my_ci_lower = ci_mod[["2.5 %"]]
-  my_ci_upper = ci_mod[["97.5 %"]]
-  
-  print(paste0('CI_lower:', my_ci_lower))
-  print(paste0('CI_upper:', my_ci_upper))
-  
-  #*************
-  # Assign the regression output in the to df (meta_pza_pnca_snps_only)
-  # you can use ('=' or '<-/->')
-  #*************
-  #my_data$logistic_logOR[my_data$mutation == i] = my_logor
- 
-  my_or -> my_data$OR[my_data$mutation == i]
-  
-  my_pval -> my_data$pvalue[my_data$mutation == i]
-  
-  my_zval -> my_data$zvalue[my_data$mutation == i]
-  
-  my_se -> my_data$logistic_se[my_data$mutation == i]
-  
-  my_ci_lower -> my_data$ci_lower[my_data$mutation == i]
-  
-  my_ci_upper -> my_data$ci_upper[my_data$mutation == i]
-  
-  #=#=#=#=#=#=#=#
-  # COMMENT: This assigns the relevant extracted output 
-  # to the df and fills NA where the mutation (row) doesn't exist 
-  # in my mutation list I am iterating over
-  #=#=#=#=#=#=#=#
-  
-}
\ No newline at end of file