added commonly used mutation format for missense muts in the gene_specific nssnp_info file
This commit is contained in:
parent
a298071309
commit
626ed3a57b
4 changed files with 194 additions and 568 deletions
|
|
@ -1,565 +0,0 @@
|
||||||
#########################################################
|
|
||||||
# TASK: To compare OR from master snps
|
|
||||||
# chisq, fisher test and logistic and adjusted logistic
|
|
||||||
#########################################################
|
|
||||||
getwd()
|
|
||||||
setwd('~/git/LSHTM_analysis/scripts')
|
|
||||||
getwd()
|
|
||||||
|
|
||||||
#install.packages("logistf")
|
|
||||||
library(logistf)
|
|
||||||
#########################################################
|
|
||||||
#%% variable assignment: input and output paths & filenames
|
|
||||||
drug = 'pyrazinamide'
|
|
||||||
gene = 'pncA'
|
|
||||||
gene_match = paste0(gene,'_p.')
|
|
||||||
cat(gene_match)
|
|
||||||
|
|
||||||
#===========
|
|
||||||
# input and output dirs
|
|
||||||
#===========
|
|
||||||
datadir = paste0('~/git/Data')
|
|
||||||
indir = paste0(datadir, '/', drug, '/', 'input')
|
|
||||||
outdir = paste0(datadir, '/', drug, '/', 'output')
|
|
||||||
|
|
||||||
#===========
|
|
||||||
# input and output files
|
|
||||||
#===========
|
|
||||||
in_filename = 'original_tanushree_data_v2.csv'
|
|
||||||
#in_filename = 'mtb_gwas_v3.csv'
|
|
||||||
infile = paste0(datadir, '/', in_filename)
|
|
||||||
cat(paste0('Reading infile1: raw snps', ' ', infile) )
|
|
||||||
|
|
||||||
# infile2: _gene associated meta snps file to extract valid snps and add calcs to.
|
|
||||||
# This is outfile3 from snps_extraction.py
|
|
||||||
in_filename_metadata = paste0(tolower(gene), '_metadata.csv')
|
|
||||||
infile_metadata = paste0(outdir, '/', in_filename_metadata)
|
|
||||||
cat(paste0('Reading infile2: gene associated metadata:', infile_metadata))
|
|
||||||
|
|
||||||
#===========
|
|
||||||
# output
|
|
||||||
#===========
|
|
||||||
out_filename = paste0(tolower(gene),'_', 'af_or.csv')
|
|
||||||
outfile = paste0(outdir, '/', out_filename)
|
|
||||||
cat(paste0('Output file with full path:', outfile))
|
|
||||||
#%% end of variable assignment for input and output files
|
|
||||||
|
|
||||||
#########################################################
|
|
||||||
# 1: Read master/raw snps stored in snps/
|
|
||||||
#####################################################
|
|
||||||
|
|
||||||
#===============
|
|
||||||
# Step 1: read raw snps (all remove entries with NA in pza column)
|
|
||||||
#===============
|
|
||||||
raw_data_all = read.csv(infile, stringsAsFactors = F)
|
|
||||||
|
|
||||||
# building cols to extract
|
|
||||||
dr_muts_col = paste0('dr_mutations_', drug)
|
|
||||||
other_muts_col = paste0('other_mutations_', drug)
|
|
||||||
|
|
||||||
cat('Extracting columns based on variables:\n'
|
|
||||||
, drug
|
|
||||||
, '\n'
|
|
||||||
, dr_muts_col
|
|
||||||
, '\n'
|
|
||||||
, other_muts_col
|
|
||||||
, '\n===============================================================')
|
|
||||||
|
|
||||||
raw_data = raw_data_all[,c("id"
|
|
||||||
, drug
|
|
||||||
, dr_muts_col
|
|
||||||
, other_muts_col)]
|
|
||||||
rm(raw_data_all)
|
|
||||||
|
|
||||||
rm(indir, in_filename, infile)
|
|
||||||
|
|
||||||
#===========
|
|
||||||
# 1a: exclude na
|
|
||||||
#===========
|
|
||||||
raw_data = raw_data[!is.na(raw_data[[drug]]),]
|
|
||||||
|
|
||||||
total_samples = length(unique(raw_data$id))
|
|
||||||
cat(paste0('Total samples without NA in', ' ', drug, 'is:', total_samples))
|
|
||||||
|
|
||||||
# sanity check: should be true
|
|
||||||
is.numeric(total_samples)
|
|
||||||
|
|
||||||
#===========
|
|
||||||
# 1b: combine the two mutation columns
|
|
||||||
#===========
|
|
||||||
#raw_data$all_mutations_pyrazinamide = paste(raw_data$dr_mutations_pyrazinamide, raw_data$other_mutations_pyrazinamide)
|
|
||||||
|
|
||||||
all_muts_colname = paste0('all_mutations_', drug)
|
|
||||||
raw_data[[all_muts_colname]] = paste(raw_data[[dr_muts_col]], raw_data[[other_muts_col]])
|
|
||||||
head(raw_data[[all_muts_colname]])
|
|
||||||
|
|
||||||
#===========
|
|
||||||
# 1c: create yet another column that contains all the mutations but in lower case
|
|
||||||
#===========
|
|
||||||
head(raw_data[[all_muts_colname]])
|
|
||||||
raw_data$all_muts_gene = tolower(raw_data[[all_muts_colname]])
|
|
||||||
head(raw_data$all_muts_gene)
|
|
||||||
|
|
||||||
# sanity checks
|
|
||||||
#table(grepl("gene_p",raw_data$all_muts_gene))
|
|
||||||
cat(paste0('converting gene match:', gene_match, ' ', 'to lowercase'))
|
|
||||||
gene_match = tolower(gene_match)
|
|
||||||
|
|
||||||
table(grepl(gene_match,raw_data$all_muts_gene))
|
|
||||||
|
|
||||||
# sanity check
|
|
||||||
if(sum(table(grepl(gene_match, raw_data$all_muts_gene))) == total_samples){
|
|
||||||
cat('PASS: Total no. of samples match')
|
|
||||||
} else{
|
|
||||||
cat('FAIL: No. of samples mismatch')
|
|
||||||
}
|
|
||||||
|
|
||||||
#########################################################
|
|
||||||
# 2: Read valid snps for which OR
|
|
||||||
# can be calculated
|
|
||||||
#########################################################
|
|
||||||
cat(paste0('Reading metadata infile:', infile_metadata))
|
|
||||||
|
|
||||||
gene_metadata = read.csv(infile_metadata
|
|
||||||
#, file.choose()
|
|
||||||
, stringsAsFactors = F
|
|
||||||
, header = T)
|
|
||||||
|
|
||||||
|
|
||||||
# clear variables
|
|
||||||
rm(in_filename_metadata, infile_metadata)
|
|
||||||
|
|
||||||
# count na in pyrazinamide column
|
|
||||||
tot_pza_na = sum(is.na(gene_metadata$pyrazinamide))
|
|
||||||
expected_rows = nrow(gene_metadata) - tot_pza_na
|
|
||||||
|
|
||||||
# drop na from the pyrazinamide colum
|
|
||||||
gene_snps_or = gene_metadata[!is.na(gene_metadata[[drug]]),]
|
|
||||||
|
|
||||||
# sanity check
|
|
||||||
if(nrow(gene_snps_or) == expected_rows){
|
|
||||||
cat('PASS: no. of rows match with expected_rows')
|
|
||||||
} else{
|
|
||||||
cat('FAIL: nrows mismatch.')
|
|
||||||
}
|
|
||||||
|
|
||||||
# extract unique snps to iterate over for AF and OR calcs
|
|
||||||
gene_snps_unique = unique(gene_snps_or$mutation)
|
|
||||||
|
|
||||||
cat(paste0('Total no. of distinct comp snps to perform OR calcs: ', length(gene_snps_unique)))
|
|
||||||
#=====================================
|
|
||||||
#OR calcs using the following 4
|
|
||||||
#1) chisq.test
|
|
||||||
#2) fisher
|
|
||||||
#3) modified chisq.test
|
|
||||||
#4) logistic
|
|
||||||
#5) adjusted logistic?
|
|
||||||
|
|
||||||
#======================================
|
|
||||||
#########################
|
|
||||||
# custom chisq function:
|
|
||||||
# To calculate OR
|
|
||||||
#########################
|
|
||||||
#i = "pnca_p.trp68gly"
|
|
||||||
#mut = grepl(i,raw_data$all_muts_gene)
|
|
||||||
#dst = raw_data[[drug]]
|
|
||||||
|
|
||||||
#x = as.numeric(mut)
|
|
||||||
#y = dst
|
|
||||||
|
|
||||||
mychisq_or = function(x,y){
|
|
||||||
tab = as.matrix(table(x,y))
|
|
||||||
a = tab[2,2]
|
|
||||||
if (a==0){ a<-0.5}
|
|
||||||
b = tab[2,1]
|
|
||||||
if (b==0){ b<-0.5}
|
|
||||||
c = tab[1,2]
|
|
||||||
if (c==0){ c<-0.5}
|
|
||||||
d = tab[1,1]
|
|
||||||
if (d==0){ d<-0.5}
|
|
||||||
(a/b)/(c/d)
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
#======================================
|
|
||||||
#==============
|
|
||||||
# TEST WITH ONE
|
|
||||||
#==============
|
|
||||||
i = "pnca_p.trp68gly"
|
|
||||||
i = "pnca_p.gln10pro"
|
|
||||||
i = "pnca_p.leu159arg"
|
|
||||||
|
|
||||||
# IV
|
|
||||||
table(grepl(i,raw_data$all_muts_gene))
|
|
||||||
mut = grepl(i,raw_data$all_muts_gene)
|
|
||||||
|
|
||||||
# DV
|
|
||||||
#dst = raw_data$pyrazinamide
|
|
||||||
dst = raw_data[[drug]] # or raw_data[,drug]
|
|
||||||
|
|
||||||
# 2X2 table
|
|
||||||
table(mut, dst)
|
|
||||||
|
|
||||||
# CV
|
|
||||||
#c = raw_data$id[mut]
|
|
||||||
c = raw_data$id[grepl(i,raw_data$all_muts_gene)]
|
|
||||||
#sid = grepl(raw_data$id[mut], raw_data$id) # warning
|
|
||||||
#argument 'pattern' has length > 1 and only the first element will be used
|
|
||||||
#grepl(raw_data$id=="ERR2512440", raw_data$id)
|
|
||||||
|
|
||||||
sid = grepl(paste(c,collapse="|"), raw_data$id)
|
|
||||||
table(sid)
|
|
||||||
|
|
||||||
# 3X2 table
|
|
||||||
table(mut, dst, sid)
|
|
||||||
|
|
||||||
#===================================================
|
|
||||||
# compare ORs from different calcs
|
|
||||||
|
|
||||||
#1) chisq
|
|
||||||
chisq.test(table(mut,dst))
|
|
||||||
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))
|
|
||||||
|
|
||||||
|
|
||||||
#2) fisher
|
|
||||||
fisher.test(table(mut, dst))
|
|
||||||
fisher.test(table(mut, dst))$p.value
|
|
||||||
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))
|
|
||||||
|
|
||||||
|
|
||||||
#3) custom chisq
|
|
||||||
or_mychisq = mychisq_or(mut,dst)
|
|
||||||
|
|
||||||
#4) logistic
|
|
||||||
summary(model<-glm(dst ~ mut, family = binomial))
|
|
||||||
|
|
||||||
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))
|
|
||||||
|
|
||||||
# extract Z of the logistic model for a given snp
|
|
||||||
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))
|
|
||||||
#logistic_ci = paste(ci_mod[["2.5 %"]], ",", 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))
|
|
||||||
|
|
||||||
# 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(paste0('Adjusted OR logistic:', or_logistic2))
|
|
||||||
pval_logistic2 = summary(model2)$coefficients[2,4]; print(paste0('Adjusted pval logistic:',pval_logistic2))
|
|
||||||
|
|
||||||
# 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)
|
|
||||||
mychisq_or(mut,dst)
|
|
||||||
})
|
|
||||||
head(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)
|
|
||||||
|
|
||||||
# allele frequency
|
|
||||||
afs = sapply(snps,function(m){
|
|
||||||
mut = grepl(m,raw_data$all_muts_gene)
|
|
||||||
mean(mut)
|
|
||||||
})
|
|
||||||
head(afs)
|
|
||||||
|
|
||||||
# logistic reg parameters: individual sapply
|
|
||||||
|
|
||||||
#--------------
|
|
||||||
## logistci or
|
|
||||||
#--------------
|
|
||||||
ors_logistic = sapply(snps,function(m){
|
|
||||||
mut = grepl(m,raw_data$all_muts_gene)
|
|
||||||
model<-glm(dst ~ mut, family = binomial)
|
|
||||||
or_logistic = exp(summary(model)$coefficients[2,1])
|
|
||||||
|
|
||||||
})
|
|
||||||
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)
|
|
||||||
pval_logistic = summary(model)$coefficients[2,4]
|
|
||||||
})
|
|
||||||
|
|
||||||
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)
|
|
||||||
logistic_se = summary(model)$coefficients[2,2]
|
|
||||||
})
|
|
||||||
|
|
||||||
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)
|
|
||||||
logistic_zval = summary(model)$coefficients[2,3]
|
|
||||||
})
|
|
||||||
|
|
||||||
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)
|
|
||||||
ci_mod = exp(confint(model))[2,]
|
|
||||||
logistic_ci_lower = ci_mod[["2.5 %"]]
|
|
||||||
})
|
|
||||||
|
|
||||||
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)
|
|
||||||
ci_mod = exp(confint(model))[2,]
|
|
||||||
logistic_ci_upper = ci_mod[["97.5 %"]]
|
|
||||||
|
|
||||||
})
|
|
||||||
|
|
||||||
head(ci_ub_logistic); head(names(ci_ub_logistic))
|
|
||||||
|
|
||||||
#--------------
|
|
||||||
# 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]
|
|
||||||
sid = grepl(paste(c,collapse="|"), raw_data$id)
|
|
||||||
model2<-glm(dst ~ mut + sid, family = binomial)
|
|
||||||
or_logistic2 = exp(summary(model2)$coefficients[2,1])
|
|
||||||
#pval_logistic2 = summary(model2)$coefficients[2,4]
|
|
||||||
})
|
|
||||||
head(logistic_ors)
|
|
||||||
|
|
||||||
#=!=!=!=!=!=!=!=!=!=!=!=!=!=!
|
|
||||||
# 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]]
|
|
||||||
|
|
||||||
# initialise an empty df
|
|
||||||
or_df = data.frame()
|
|
||||||
|
|
||||||
x = sapply(snps,function(m){
|
|
||||||
|
|
||||||
#df = data.frame()
|
|
||||||
|
|
||||||
mut = grepl(m,raw_data$all_muts_gene)
|
|
||||||
model<-glm(dst ~ mut, family = binomial)
|
|
||||||
|
|
||||||
# allele frequency
|
|
||||||
afs = mean(mut)
|
|
||||||
|
|
||||||
# logistic model
|
|
||||||
beta_logistic = summary(model)$coefficients[2,1]
|
|
||||||
|
|
||||||
or_logistic = exp(summary(model)$coefficients[2,1])
|
|
||||||
print(paste0('logistic OR:', or_logistic))
|
|
||||||
|
|
||||||
pval_logistic = summary(model)$coefficients[2,4]
|
|
||||||
print(paste0('logistic pval:', pval_logistic))
|
|
||||||
|
|
||||||
se_logistic = summary(model)$coefficients[2,2]
|
|
||||||
zval_logistic = summary(model)$coefficients[2,3]
|
|
||||||
ci_mod = exp(confint(model))[2,]
|
|
||||||
ci_lower_logistic = ci_mod[["2.5 %"]]
|
|
||||||
ci_upper_logistic = ci_mod[["97.5 %"]]
|
|
||||||
|
|
||||||
# custom_chisq and fisher: OR p-value and CI
|
|
||||||
or_mychisq = mychisq_or(dst, mut)
|
|
||||||
|
|
||||||
or_fisher = fisher.test(dst, mut)$estimate
|
|
||||||
or_fisher = or_fisher[[1]]
|
|
||||||
|
|
||||||
pval_fisher = fisher.test(dst, mut)$p.value
|
|
||||||
|
|
||||||
ci_lower_fisher = fisher.test(dst, mut)$conf.int[1]
|
|
||||||
ci_upper_fisher = fisher.test(dst, mut)$conf.int[2]
|
|
||||||
|
|
||||||
# chi sq estimates
|
|
||||||
estimate_chisq = chisq.test(dst, mut)$statistic; estimate_chisq
|
|
||||||
est_chisq = estimate_chisq[[1]]; print(est_chisq)
|
|
||||||
|
|
||||||
pval_chisq = chisq.test(dst, mut)$p.value
|
|
||||||
|
|
||||||
# build a row to append to df
|
|
||||||
row = data.frame(mutation = m
|
|
||||||
, af = afs
|
|
||||||
, beta_logistic = beta_logistic
|
|
||||||
, or_logistic = or_logistic
|
|
||||||
, pval_logistic = pval_logistic
|
|
||||||
, se_logistic = se_logistic
|
|
||||||
, zval_logistic = zval_logistic
|
|
||||||
, ci_low_logistic = ci_lower_logistic
|
|
||||||
, ci_hi_logistic = ci_upper_logistic
|
|
||||||
, or_mychisq = or_mychisq
|
|
||||||
, or_fisher = or_fisher
|
|
||||||
, pval_fisher = pval_fisher
|
|
||||||
, ci_low_fisher= ci_lower_fisher
|
|
||||||
, ci_hi_fisher = ci_upper_fisher
|
|
||||||
, est_chisq = est_chisq
|
|
||||||
, pval_chisq = pval_chisq
|
|
||||||
)
|
|
||||||
#print(row)
|
|
||||||
|
|
||||||
or_df <<- rbind(or_df, row)
|
|
||||||
|
|
||||||
})
|
|
||||||
|
|
||||||
write.csv(or_df, 'test_ors.csv')
|
|
||||||
|
|
||||||
#=================================
|
|
||||||
# testing logistic or with maxit, etc.
|
|
||||||
|
|
||||||
print(paste0('subset to iterate over;', snps))
|
|
||||||
|
|
||||||
# start loop
|
|
||||||
perfectSeparation <- function(w) {
|
|
||||||
if(grepl("fitted probabilities numerically 0 or 1 occurred",
|
|
||||||
as.character(w))) {ww <<- ww+1}
|
|
||||||
}
|
|
||||||
|
|
||||||
for (i in snps){
|
|
||||||
|
|
||||||
print(i)
|
|
||||||
|
|
||||||
# IV
|
|
||||||
#mut<-as.numeric(grepl(i,raw_data$all_muts_gene))
|
|
||||||
mut = grepl(i,raw_data$all_muts_gene)
|
|
||||||
table(mut)
|
|
||||||
|
|
||||||
# DV
|
|
||||||
#dst<-as.numeric(raw_data[[drug]])
|
|
||||||
dst = raw_data[[drug]]
|
|
||||||
|
|
||||||
# table
|
|
||||||
print(table(dst, mut))
|
|
||||||
|
|
||||||
#=====================
|
|
||||||
# logistic regression, glm.control(maxit = n)
|
|
||||||
#https://stats.stackexchange.com/questions/11109/how-to-deal-with-perfect-separation-in-logistic-regression
|
|
||||||
#=====================
|
|
||||||
#n = 1
|
|
||||||
summary(model<-glm(dst ~ mut
|
|
||||||
, family = binomial
|
|
||||||
#, control = glm.control(maxit = 1)
|
|
||||||
#, options(warn = 1)
|
|
||||||
))
|
|
||||||
#, warning = perfectSeparation))
|
|
||||||
or_logistic = exp(summary(model)$coefficients[2,1]); print(or_logistic)
|
|
||||||
pval_logistic_maxit = summary(model)$coefficients[2,4]; print(pval_logistic_maxit)
|
|
||||||
logistic_se = summary(model)$coefficients[2,2]
|
|
||||||
logistic_zval = summary(model)$coefficients[2,3]
|
|
||||||
ci_mod = exp(confint(model))[2,]
|
|
||||||
logistic_ci_lower = ci_mod[["2.5 %"]]
|
|
||||||
logistic_ci_upper = ci_mod[["97.5 %"]]
|
|
||||||
|
|
||||||
#=====================
|
|
||||||
# fishers test
|
|
||||||
#=====================
|
|
||||||
#attributes(fisher.test(table(dst, mut)))
|
|
||||||
or_fisher = fisher.test(table(dst, mut))$estimate
|
|
||||||
or_fisher = or_fisher[[1]]; or_fisher
|
|
||||||
|
|
||||||
pval_fisher = fisher.test(table(dst, mut))$p.value ; pval_fisher
|
|
||||||
|
|
||||||
#=====================
|
|
||||||
# custom chi square
|
|
||||||
#=====================
|
|
||||||
or_mychisq = mychisq_or(mut,dst)
|
|
||||||
|
|
||||||
#=====================
|
|
||||||
# chi square
|
|
||||||
#=====================
|
|
||||||
#chisq.test(y = dst, x = mut)
|
|
||||||
#attributes(chisq.test(table(dst, mut)))
|
|
||||||
est_chisq = chisq.test(table(dst, mut))$statistic
|
|
||||||
est_chisq = est_chisq[[1]]; est_chisq
|
|
||||||
|
|
||||||
pval_chisq = chisq.test(table(dst, mut))$p.value; pval_chisq
|
|
||||||
|
|
||||||
# all output
|
|
||||||
writeLines(c(paste0("mutation:", i)
|
|
||||||
, paste0("=========================")
|
|
||||||
, 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
|
|
||||||
#=====================
|
|
||||||
#attributes(fisher.test(table(dst, mut)))
|
|
||||||
or_fisher = fisher.test(table(dst, mut))$estimate
|
|
||||||
or_fisher = or_fisher[[1]]; or_fisher
|
|
||||||
|
|
||||||
pval_fisher = fisher.test(table(dst, mut))$p.value ; pval_fisher
|
|
||||||
|
|
||||||
|
|
||||||
# https://stats.stackexchange.com/questions/259635/what-is-the-difference-using-a-fishers-exact-test-vs-a-logistic-regression-for
|
|
||||||
#mymethod = "minlike" # default
|
|
||||||
mymethod = "central"
|
|
||||||
#mymethod = "blaker"
|
|
||||||
|
|
||||||
exact2x2(table(dst, mut),tsmethod=mymethod)
|
|
||||||
|
|
||||||
mcnemar.exact(x,y=NULL, conf.level=.95)
|
|
||||||
|
|
||||||
#=====================================================================
|
|
||||||
152
scripts/nssnp_info_format.py
Executable file
152
scripts/nssnp_info_format.py
Executable file
|
|
@ -0,0 +1,152 @@
|
||||||
|
#!/usr/bin/env python3
|
||||||
|
# -*- coding: utf-8 -*-
|
||||||
|
"""
|
||||||
|
Created on Wed Jun 10 11:13:49 2020
|
||||||
|
|
||||||
|
@author: tanu
|
||||||
|
"""
|
||||||
|
#==============================================================================
|
||||||
|
# TASK
|
||||||
|
# To format snp_fino.txt file that has already been processed in bash
|
||||||
|
# to include mcsm style muts and gwas style muts. The idea is that the info file
|
||||||
|
# will contain all possible mutation format style to make it easy to populate
|
||||||
|
# and link other files with this sort of meta data. For example: or file
|
||||||
|
#=======================================================================
|
||||||
|
|
||||||
|
# FIXME : add bash info here as well
|
||||||
|
|
||||||
|
#%% useful links
|
||||||
|
#https://chrisalbon.com/python/data_wrangling/pandas_join_merge_dataframe/
|
||||||
|
#https://kanoki.org/2019/11/12/how-to-use-regex-in-pandas/
|
||||||
|
#https://stackoverflow.com/questions/40348541/pandas-diff-with-string
|
||||||
|
#=======================================================================
|
||||||
|
#%% specify dirs
|
||||||
|
import os, sys
|
||||||
|
import pandas as pd
|
||||||
|
import numpy as np
|
||||||
|
import re
|
||||||
|
import argparse
|
||||||
|
|
||||||
|
homedir = os.path.expanduser('~')
|
||||||
|
os.chdir(homedir + '/git/LSHTM_analysis/scripts')
|
||||||
|
|
||||||
|
#from reference_dict import my_aa_dict
|
||||||
|
#from reference_dict import low_3letter_dict # equivalent of my_aa_dict
|
||||||
|
from reference_dict import oneletter_aa_dict
|
||||||
|
#=======================================================================
|
||||||
|
#%% command line args
|
||||||
|
arg_parser = argparse.ArgumentParser()
|
||||||
|
|
||||||
|
arg_parser.add_argument('-d', '--drug', help = 'drug name', default = 'pyrazinamide')
|
||||||
|
arg_parser.add_argument('-g', '--gene', help = 'gene name (case sensitive)', default = 'pncA') # case sensitive
|
||||||
|
|
||||||
|
args = arg_parser.parse_args()
|
||||||
|
#=======================================================================
|
||||||
|
#%% variables
|
||||||
|
#gene = 'pncA'
|
||||||
|
#drug = 'pyrazinamide'
|
||||||
|
#gene_match = gene +'_p.'
|
||||||
|
|
||||||
|
# cmd variables
|
||||||
|
gene = args.gene
|
||||||
|
drug = args.drug
|
||||||
|
gene_match = gene +'_p.'
|
||||||
|
|
||||||
|
#=======================================================================
|
||||||
|
#%% input and output dirs and files
|
||||||
|
#=======
|
||||||
|
# data dir
|
||||||
|
#=======
|
||||||
|
datadir = homedir + '/' + 'git/Data'
|
||||||
|
indir = datadir + '/' + drug + '/input'
|
||||||
|
outdir = datadir + '/' + drug + '/output'
|
||||||
|
|
||||||
|
#=======
|
||||||
|
# input
|
||||||
|
#=======
|
||||||
|
|
||||||
|
gene_info_filename = 'ns'+ gene.lower()+ '_snp_info1.txt'
|
||||||
|
gene_info = indir + '/' + gene_info_filename
|
||||||
|
print('gene info file: ', gene_info
|
||||||
|
, '\n============================================================')
|
||||||
|
|
||||||
|
#=======
|
||||||
|
# output
|
||||||
|
#=======
|
||||||
|
gene_snp_info_filename = 'ns' + gene.lower() + '_snp_info.csv' # other one is called AFandOR
|
||||||
|
outfile_snp_info = indir + '/' + gene_snp_info_filename
|
||||||
|
print('Output file: ', outfile_snp_info
|
||||||
|
, '\n============================================================')
|
||||||
|
|
||||||
|
#%% read files: preformatted using bash
|
||||||
|
info_df2 = pd.read_csv(gene_info, sep = '\t', header = 0) #303, 10
|
||||||
|
|
||||||
|
#%% Split into three cols with 1-letter aa_code & combine to get mutationinformation column
|
||||||
|
# check mutation format in exisiting df
|
||||||
|
info_df2.head()
|
||||||
|
info_df2['mut_info'].head()
|
||||||
|
|
||||||
|
print('Creating column: mutationinformation')
|
||||||
|
info_df2_ncols = len(info_df2.columns)
|
||||||
|
|
||||||
|
info_df2['wild_type'] = info_df2['mut_info'].str.extract('(\w{1})>')
|
||||||
|
info_df2['position'] = info_df2['mut_info'].str.extract('(\d+)')
|
||||||
|
info_df2['mutant_type'] = info_df2['mut_info'].str.extract('>\d+(\w{1})')
|
||||||
|
|
||||||
|
info_df2['mutationinformation'] = info_df2['wild_type'] + info_df2['position'] + info_df2['mutant_type']
|
||||||
|
|
||||||
|
# sanity check
|
||||||
|
ncols_add = 4 # Beware hardcoded
|
||||||
|
if len(info_df2.columns) == info_df2_ncols + ncols_add:
|
||||||
|
print('PASS: Succesfully extracted and added mutationinformation (mcsm style) as below\n'
|
||||||
|
, info_df2['mutationinformation'].head()
|
||||||
|
, '\n=====================================================================')
|
||||||
|
else:
|
||||||
|
print('FAIL: No. of cols mismatch'
|
||||||
|
,'\noriginal length:', info_df2_ncols
|
||||||
|
, '\nExpected no. of cols:', info_df2_ncols + ncols_add
|
||||||
|
, '\nGot no. of cols:', len(info_df2.columns))
|
||||||
|
sys.exit()
|
||||||
|
|
||||||
|
# update ncols
|
||||||
|
info_df2_ncols = len(info_df2.columns)
|
||||||
|
|
||||||
|
#%% Creating column 'mutation' which as mutation of the format;
|
||||||
|
# <gene_match>.lower()<abc>1<cde>: pnca_p.trp68gly
|
||||||
|
# match the 'one_letter_code' value to get the dict key (three-letter code)
|
||||||
|
|
||||||
|
print('Creating column: mutation')
|
||||||
|
|
||||||
|
# dict to use: oneletter_aa_dict
|
||||||
|
lookup_dict = dict()
|
||||||
|
for k1, v1 in oneletter_aa_dict.items():
|
||||||
|
lookup_dict[k1] = v1['three_letter_code_lower']
|
||||||
|
for k2, v2 in lookup_dict.items():
|
||||||
|
info_df2['wt_3let'] = info_df2['wild_type'].squeeze().map(lookup_dict)
|
||||||
|
info_df2['mt_3let'] = info_df2['mutant_type'].squeeze().map(lookup_dict)
|
||||||
|
|
||||||
|
# create column mutation
|
||||||
|
info_df2['mutation'] = info_df2['wt_3let'] + info_df2['position'] + info_df2['mt_3let']
|
||||||
|
|
||||||
|
# add prefix: gene_match to each value in column
|
||||||
|
info_df2['mutation'] = gene_match.lower() + info_df2['mutation'].astype(str)
|
||||||
|
|
||||||
|
# sanity check
|
||||||
|
ncols_add = 3 # Beware hardcoded
|
||||||
|
if len(info_df2.columns) == info_df2_ncols + ncols_add:
|
||||||
|
print('PASS: Succesfully created column mutation as below\n'
|
||||||
|
, info_df2['mutation'].head()
|
||||||
|
, '\n=====================================================================')
|
||||||
|
else:
|
||||||
|
print('FAIL: No. of cols mismatch\noriginal length:', info_df2_ncols
|
||||||
|
, '\nExpected no. of cols:', info_df2_ncols + ncols_add
|
||||||
|
, '\nGot no. of cols:', len(info_df2.columns))
|
||||||
|
sys.exit()
|
||||||
|
|
||||||
|
#%% write file
|
||||||
|
print('\n====================================================================='
|
||||||
|
, '\nWriting output file:\n', outfile_snp_info
|
||||||
|
, '\nNo.of rows:', len(info_df2)
|
||||||
|
, '\nNo. of cols:', len(info_df2.columns))
|
||||||
|
info_df2.to_csv(outfile_snp_info, index = False)
|
||||||
|
|
||||||
|
|
@ -221,9 +221,9 @@ else:
|
||||||
print('FAIL: Second cross check unsuccessfull. Debug please!')
|
print('FAIL: Second cross check unsuccessfull. Debug please!')
|
||||||
sys.exit()
|
sys.exit()
|
||||||
|
|
||||||
|
#%% extract mut info into three cols
|
||||||
orig_len = len(dfm2_mis.columns)
|
orig_len = len(dfm2_mis.columns)
|
||||||
|
|
||||||
#%% extract mut info into three cols
|
|
||||||
dfm2_mis['wild_type'] = dfm2_mis['mut_info'].str.extract('(\w{1})>')
|
dfm2_mis['wild_type'] = dfm2_mis['mut_info'].str.extract('(\w{1})>')
|
||||||
dfm2_mis['position'] = dfm2_mis['mut_info'].str.extract('(\d+)')
|
dfm2_mis['position'] = dfm2_mis['mut_info'].str.extract('(\d+)')
|
||||||
dfm2_mis['mutant_type'] = dfm2_mis['mut_info'].str.extract('>\d+(\w{1})')
|
dfm2_mis['mutant_type'] = dfm2_mis['mut_info'].str.extract('>\d+(\w{1})')
|
||||||
|
|
|
||||||
|
|
@ -57,12 +57,17 @@ print('Input filename:', in_filename
|
||||||
#%% end of variable assignment for input and output files
|
#%% end of variable assignment for input and output files
|
||||||
#=======================================================================
|
#=======================================================================
|
||||||
#%% Read input file
|
#%% Read input file
|
||||||
my_aa = pd.read_csv(infile) #20, 6
|
aa_table = pd.read_csv(infile) #20, 6
|
||||||
|
|
||||||
|
#------------------------
|
||||||
|
#1) 3-letter (lower) code as key
|
||||||
|
#-------------------------
|
||||||
# assign the one_letter code as the row names so that it is easier to create
|
# assign the one_letter code as the row names so that it is easier to create
|
||||||
# a dict of dicts using index
|
# a dict of dicts using index
|
||||||
#my_aa = pd.read_csv('aa_codes.csv', index_col = 0) #20, 6 #a way to it since it is the first column
|
#my_aa = pd.read_csv('aa_codes.csv', index_col = 0) #20, 6 #a way to it since it is the first column
|
||||||
my_aa = my_aa.set_index('three_letter_code_lower') #20, 5
|
my_aa = aa_table.set_index('three_letter_code_lower') #20, 5
|
||||||
|
my_aa.columns
|
||||||
|
my_aa.index
|
||||||
|
|
||||||
#==================
|
#==================
|
||||||
# convert file
|
# convert file
|
||||||
|
|
@ -75,6 +80,40 @@ my_aa = my_aa.set_index('three_letter_code_lower') #20, 5
|
||||||
my_aa_dict = my_aa.to_dict('index') #20, with 5 subkeys
|
my_aa_dict = my_aa.to_dict('index') #20, with 5 subkeys
|
||||||
print('Printing my_aa_dict:', my_aa_dict.keys())
|
print('Printing my_aa_dict:', my_aa_dict.keys())
|
||||||
|
|
||||||
|
#FIXME : use the below in all code
|
||||||
|
low_3letter_dict = my_aa.to_dict('index') #20, with 5 subkeys
|
||||||
|
print('Printing lower-case 3 letter aa dict:',low_3letter_dict.keys())
|
||||||
|
|
||||||
|
#------------------------
|
||||||
|
#2) 1-letter code as key
|
||||||
|
#-------------------------
|
||||||
|
aa_1let = aa_table.set_index('one_letter_code') #20, 5
|
||||||
|
aa_1let.columns
|
||||||
|
aa_1let.index
|
||||||
|
|
||||||
|
oneletter_aa_dict = aa_1let.to_dict('index') #20, with 5 subkeys
|
||||||
|
print('Printing one letter aa dict:', oneletter_aa_dict.keys())
|
||||||
|
|
||||||
|
#------------------------
|
||||||
|
#3) amino acid name as key
|
||||||
|
#-------------------------
|
||||||
|
aa_name = aa_table.set_index('amino_acid_name') #20, 5
|
||||||
|
aa_name.columns
|
||||||
|
aa_name.index
|
||||||
|
|
||||||
|
aa_name_dict = aa_name.to_dict('index') #20, with 5 subkeys
|
||||||
|
print('Printing amino acid names aa dict:', aa_name_dict.keys())
|
||||||
|
|
||||||
|
#------------------------
|
||||||
|
#3) 3 letter uppercase as key
|
||||||
|
#-------------------------
|
||||||
|
aa_up3let = aa_table.set_index('three_letter_code_upper') #20, 5
|
||||||
|
aa_up3let.columns
|
||||||
|
aa_up3let.index
|
||||||
|
|
||||||
|
up_3letter_aa_dict = aa_up3let.to_dict('index') #20, with 5 subkeys
|
||||||
|
print('Printing upper case 3 letter aa dict:', up_3letter_aa_dict.keys())
|
||||||
|
|
||||||
#================================================
|
#================================================
|
||||||
# dict of aa with their corresponding properties
|
# dict of aa with their corresponding properties
|
||||||
# This is defined twice
|
# This is defined twice
|
||||||
|
|
|
||||||
Loading…
Add table
Add a link
Reference in a new issue