tanu/LSHTM_analysis
1
0
Fork 0
Code Issues Pull requests Releases Wiki Activity

tody scracth script for various OR calcs

Browse source
This commit is contained in:
Tanushree Tunstall 2020-06-23 11:57:51 +01:00
parent 18998092f4
commit 1e43ca8136
1 changed files with 104 additions and 184 deletions
Download patch file Download diff file Expand all files Collapse all files

288
scripts/af_or_calcs_scratch.R
View file

@ -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
#=#=#=#=#=#=#=#
}