241 lines
6.5 KiB
R
241 lines
6.5 KiB
R
getwd()
|
|
setwd("/git/github/git/LSHTM_analysis/meta_data_analysis")
|
|
getwd()
|
|
|
|
#===============
|
|
# Step 1: read GWAS raw data stored in Data_original/
|
|
#===============
|
|
infile = read.csv(file.choose(), stringsAsFactors = F)
|
|
|
|
raw_data = infile[,c("id"
|
|
, "pyrazinamide"
|
|
, "dr_mutations_pyrazinamide"
|
|
, "other_mutations_pyrazinamide")]
|
|
|
|
#####
|
|
# 1a: exclude na
|
|
#####
|
|
raw_data = raw_data[!is.na(raw_data$pyrazinamide),]
|
|
|
|
total_samples = length(unique(raw_data$id))
|
|
print(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)
|
|
head(raw_data$all_mutations_pyrazinamide)
|
|
|
|
#####
|
|
# 1c: create yet another column that contains all the mutations but in lower case
|
|
#####
|
|
raw_data$all_muts_pnca = tolower(raw_data$all_mutations_pyrazinamide)
|
|
|
|
# sanity checks
|
|
table(grepl("pnca_p",raw_data$all_muts_pnca))
|
|
|
|
# sanity check: should be TRUE
|
|
sum(table(grepl("pnca_p",raw_data$all_muts_pnca))) == total_samples
|
|
|
|
# set up variables: can be used for logistic regression as well
|
|
i = "pnca_p.ala134gly" # has a NA, should NOT exist
|
|
table(grepl(i,raw_data$all_muts_pnca))
|
|
|
|
i = "pnca_p.trp68gly"
|
|
table(grepl(i,raw_data$all_muts_pnca))
|
|
|
|
mut = grepl(i,raw_data$all_muts_pnca)
|
|
dst = raw_data$pyrazinamide
|
|
table(mut, dst)
|
|
|
|
#chisq.test(table(mut,dst))
|
|
#fisher.test(table(mut, dst))
|
|
#table(mut)
|
|
|
|
###### read list of muts to calculate OR for (fname3 from pnca_data_extraction.py)
|
|
pnca_snps_or = read.csv(file.choose()
|
|
, stringsAsFactors = F
|
|
, header = T)
|
|
|
|
# extract unique snps to iterate over for AF and OR calcs
|
|
# total no of unique snps
|
|
# AF and OR calculations
|
|
|
|
pnca_snps_unique = unique(pnca_snps_or$mutation)
|
|
|
|
# Define OR function
|
|
x = as.numeric(mut)
|
|
y = dst
|
|
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)
|
|
|
|
}
|
|
|
|
dst = raw_data$pyrazinamide
|
|
ors = sapply(pnca_snps_unique,function(m){
|
|
mut = grepl(m,raw_data$all_muts_pnca)
|
|
or(mut,dst)
|
|
})
|
|
|
|
ors
|
|
|
|
pvals = sapply(pnca_snps_unique,function(m){
|
|
mut = grepl(m,raw_data$all_muts_pnca)
|
|
fisher.test(mut,dst)$p.value
|
|
})
|
|
|
|
pvals
|
|
|
|
afs = sapply(pnca_snps_unique,function(m){
|
|
mut = grepl(m,raw_data$all_muts_pnca)
|
|
mean(mut)
|
|
})
|
|
|
|
afs
|
|
|
|
# check ..hmmm
|
|
afs['pnca_p.trp68gly']
|
|
afs['pnca_p.gln10pro']
|
|
afs['pnca_p.leu4ser']
|
|
|
|
#plot(density(log(ors)))
|
|
#plot(-log10(pvals))
|
|
#hist(log(ors)
|
|
# ,breaks = 100
|
|
# )
|
|
|
|
# subset df cols to add to the calc param df
|
|
pnca_snps_cols = pnca_snps_or[c('mutation_info', 'mutation', 'Mutationinformation')]
|
|
pnca_snps_cols = pnca_snps_cols[!duplicated(pnca_snps_cols$mutation),]
|
|
|
|
rownames(pnca_snps_cols) = pnca_snps_cols$mutation
|
|
head(rownames(pnca_snps_cols))
|
|
#snps_with_AF_and_OR
|
|
|
|
# combine
|
|
comb_AF_and_OR = data.frame(ors, pvals, afs)
|
|
head(rownames(comb_AF_and_OR))
|
|
|
|
# sanity checks: should be the same
|
|
dim(comb_AF_and_OR); dim(pnca_snps_cols)
|
|
table(rownames(comb_AF_and_OR)%in%rownames(pnca_snps_cols))
|
|
|
|
table(rownames(pnca_snps_cols)%in%rownames(comb_AF_and_OR))
|
|
|
|
# merge the above two df whose dim you checked
|
|
snps_with_AF_and_OR = merge(comb_AF_and_OR, pnca_snps_cols
|
|
, by = "row.names"
|
|
# , all.x = T
|
|
)
|
|
|
|
#rm(pnca_snps_cols, pnca_snps_or, raw_data)
|
|
|
|
#===============
|
|
# Step 3: Read data file where you will add the calculated OR
|
|
# Note: this is the big file with one-many relationship between snps and lineages
|
|
# i.e fname4 from 'pnca_extraction.py'
|
|
#===============
|
|
my_data = read.csv(file.choose()
|
|
, row.names = 1
|
|
, stringsAsFactors = FALSE)
|
|
|
|
head(my_data)
|
|
length(unique(my_data$id))
|
|
|
|
# check if first col is 'id': should be TRUE
|
|
colnames(my_data)[1] == 'id'
|
|
|
|
# sanity check
|
|
head(my_data$mutation)
|
|
|
|
# FILES TO MERGE:
|
|
# comb_AF_and_OR: file containing OR
|
|
# my_data = big meta data file
|
|
# linking column: mutation
|
|
|
|
head(my_data)
|
|
merged_df = merge(my_data # big file
|
|
, snps_with_AF_and_OR # small (afor file)
|
|
, by = "mutation"
|
|
, all.x = T) # because you want all the entries of the meta data
|
|
|
|
# sanity checks: should be True
|
|
# FIXME: I have checked this manually, but make it so it is a pass or a fail!
|
|
comb_AF_and_OR[rownames(comb_AF_and_OR) == "pnca_p.gln10pro",]$ors
|
|
merged_df[merged_df$Mutationinformation.x == "Q10P",]$ors
|
|
|
|
merged_df[merged_df$Mutationinformation.x == "Q10P",]
|
|
|
|
# sanity check: very important!
|
|
colnames(merged_df)
|
|
|
|
table(merged_df$mutation_info.x == merged_df$mutation_info.y)
|
|
|
|
#FIXME: what happened to other 7 and FALSE
|
|
table(merged_df$Mutationinformation.x == merged_df$Mutationinformation.y)
|
|
|
|
# problem
|
|
identical(merged_df$Mutationinformation.x, merged_df$Mutationinformation.y)
|
|
|
|
#merged_df[merged_df$Mutationinformation.x != merged_df$Mutationinformation.y,]
|
|
|
|
#throw away the y because that is a smaller df
|
|
d1 = which(colnames(merged_df) == "mutation_info.y") #21
|
|
d2 = which(colnames(merged_df) == "Mutationinformation.y") #22
|
|
|
|
merged_df2 = merged_df[-c(d1, d2)] #3093 20
|
|
colnames(merged_df2)
|
|
|
|
# rename cols
|
|
colnames(merged_df2)[colnames(merged_df2)== "mutation_info.x"] <- "mutation_info"
|
|
colnames(merged_df2)[colnames(merged_df2)== "Mutationinformation.x"] <- "Mutationinformation"
|
|
|
|
colnames(merged_df2)
|
|
|
|
# should be 0
|
|
sum(is.na(merged_df2$Mutationinformation))
|
|
|
|
# count na in each column
|
|
na_count = sapply(merged_df2, function(y) sum(length(which(is.na(y))))); na_count
|
|
# only some or and Af should be NA
|
|
#Row.names ors pvals afs
|
|
#81 81 81 81
|
|
|
|
|
|
colnames(merged_df2)[colnames(merged_df2)== "ors"] <- "OR"
|
|
colnames(merged_df2)[colnames(merged_df2)== "afs"] <- "AF"
|
|
colnames(merged_df2)[colnames(merged_df2)== "pvals"] <- "pvalue"
|
|
|
|
colnames(merged_df2)
|
|
|
|
# add log OR and neglog pvalue
|
|
merged_df2$logor = log(merged_df2$OR)
|
|
is.numeric(merged_df2$logor)
|
|
|
|
merged_df2$neglog10pvalue = -log10(merged_df2$pvalue)
|
|
is.numeric(merged_df2$neglog10pvalue)
|
|
|
|
# write file out
|
|
#write.csv(merged_df, "../Data/meta_data_with_AFandOR_JP_TT.csv")
|
|
|
|
# define output variables
|
|
drug = 'pyrazinamide'
|
|
out_dir = paste0("../mcsm_analysis/",drug,"/Data/")
|
|
outFile = "meta_data_with_AFandOR.csv"
|
|
output_filename = paste0(outdir, outFile)
|
|
|
|
write.csv(merged_df2, output_filename
|
|
, row.names = F)
|