suppressMessages(silent <- lapply(
c("biomaRt", "tidyverse","vroom", "ggrastr", "hudson"),
library, character.only=T))
table = function (..., useNA = 'always') base::table(..., useNA = useNA)p="Veg5yr1"
dir="/Users/mike/Documents/Research/Vegetarian/regenie-GWAS/sumstats-raw/noBMI/"
outdir="Interactive.GWAS"
chr=list()
for (c in 1:22){
chr[[c]]=vroom(paste(
dir, "/Veg5yr1_GWAS_chr",c,"_",p,".regenie", sep=""), show_col_types = FALSE,progress = FALSE)
}
file=do.call(rbind,chr)
dirbmi="/Users/mike/Documents/Research/Vegetarian/regenie-GWAS/sumstats-raw/adjBMI"
bmi=list()
for (c in 1:22){
bmi[[c]]=vroom(paste(
dirbmi, "/Veg5yr1_GWAS_chr",c,"_",p,".regenie", sep=""), show_col_types = FALSE,progress = FALSE)
}
filebmi=do.call(rbind,bmi)
# interactive GWAS
file$P=10^(-file$LOG10P)
file=file%>%select(-LOG10P)
colnames(file)=c("CHR", "POS", "SNP",
"ALLELE0","ALLELE1","A1FREQ","INFO" ,"N",
"TEST","BETA","SE","CHISQ","EXTRA","pvalue")
file$Link=ifelse(file$pvalue<5e-05,
paste0("https://www.ncbi.nlm.nih.gov/snp/",
file$SNP),
NA)
file$Hover= ifelse(file$pvalue<5e-05,
paste0("SNP: ", file$SNP,
"\nPOS: ", file$CHR, ":", file$POS,
"\nAllele: ", file$ALLELE1,
"\np-value: ", formatC(file$pvalue, format = "e", digits = 3),
"\nBeta (SE): ", paste(signif(filebmi$BETA, digits = 4), "(", signif(filebmi$SE, digits = 4),")"),
"\nN: ", file$N),
NA)
filebmi$P=10^(-filebmi$LOG10P)
filebmi=filebmi%>%select(-LOG10P)
colnames(filebmi)=c("CHR", "POS", "SNP",
"ALLELE0","ALLELE1","A1FREQ","INFO" ,"N",
"TEST","BETA","SE","CHISQ","EXTRA","pvalue")
filebmi$Link=ifelse(filebmi$pvalue<5e-05,
paste0("https://www.ncbi.nlm.nih.gov/snp/",
filebmi$SNP),
NA)
filebmi$Hover= ifelse(filebmi$pvalue<5e-05,
paste0("SNP: ", filebmi$SNP,
"\nPOS: ", filebmi$CHR, ":", filebmi$POS,
"\nAllele: ", filebmi$ALLELE1,
"\np-value: ", formatC(filebmi$pvalue, format = "e", digits = 3),
"\nBeta (SE): ", paste(signif(filebmi$BETA, digits = 4), "(", signif(filebmi$SE, digits = 4),")"),
"\nN: ", filebmi$N),
NA)
igmirror(top=file, bottom=filebmi,
tline=5e-08,
bline=5e-08,
toptitle="GWAS of Vegetarianism (no BMI in model)",
bottomtitle = "GWAS of Vegetarianism (BMI adjusted)",
highlight_p = c(5e-05, 5e-05),
highlighter="green",
file=paste0(outdir, "/PlotVegGWAS.html")
)dir="/Users/mike/Documents/Research/Vegetarian/regenie-GWAS/MAGMA/sumstats"
dirncbi="/Users/mike/Documents/Research/Vegetarian/MAGMA5/sumstats/noBMI"
outdir="Interactive.GWAS.gene"
dir.create(outdir)
#Prep gene annotation file
geneloc=as_tibble(read.table(paste0(dirncbi, "/NCBI37.3.gene.loc")))
colnames(geneloc)=c("GENE","CHR","START","STOP", "STRAND", "GENE_NAME")
geneloc=geneloc%>%filter(CHR!="X" & CHR!="Y")
geneloc$CHR=as.integer(geneloc$CHR)
geneloc=geneloc%>%select(GENE, GENE_NAME)
file=read_table(paste0(dir, "/2.geneanalysis.genes.out"))
filebmi=read_table(paste0(dir, "/2.geneanalysis.BMI.genes.out"))
#only use ID and name as start and stop will be different because of the +2kb and -1kb window added during analysis
file=left_join(file,geneloc)
file=file%>%select(GENE_NAME, everything())
filebmi=left_join(filebmi,geneloc)
filebmi=filebmi%>%select(GENE_NAME, everything())
#Set colnames to hudson format
colnames(file)=c("SNP","GENE_ID_NUMBER", "CHR", "POS",
"STOP","NSNPS","NPARAM",
"N","ZSTAT","pvalue",
"P_SNPWISE_MEAN", "P_SNPWISE_TOP1")
colnames(filebmi)=c("SNP","GENE_ID_NUMBER", "CHR", "POS",
"STOP","NSNPS","NPARAM",
"N","ZSTAT","pvalue",
"P_SNPWISE_MEAN", "P_SNPWISE_TOP1")
#Add Link col
file$Link=ifelse(file$pvalue<1e-3,paste0("https://www.ncbi.nlm.nih.gov/gene/?term=", file$SNP),NA)
filebmi$Link=ifelse(file$pvalue<1e-3,paste0("https://www.ncbi.nlm.nih.gov/gene/?term=", filebmi$SNP),NA)
#Add Hover col
file$Hover= ifelse(file$pvalue<1e-3,paste0("Gene: ", file$SNP,
"\nStart-Stop: ", file$CHR, ":", file$POS, "-", file$STOP,
"\nNSNPS: ", file$NSNPS,
"\nNPARAM: ", file$NPARAM,
"\nP_MULTI: ", formatC(file$pvalue, format = "e", digits = 3),
"\nP_SNPWISE_MEAN: ", formatC(file$P_SNPWISE_MEAN, format = "e", digits = 3),
"\nP_SNPWISE_TOP1: ", formatC(file$P_SNPWISE_TOP1, format = "e", digits = 3),
"\nN: ", file$N), NA)
filebmi$Hover= ifelse(filebmi$pvalue<1e-3,paste0("Gene: ", filebmi$SNP,
"\nStart-Stop: ", filebmi$CHR, ":", filebmi$POS, "-", filebmi$STOP,
"\nNSNPS: ", filebmi$NSNPS,
"\nNPARAM: ", filebmi$NPARAM,
"\nP_MULTI: ", formatC(filebmi$pvalue, format = "e", digits = 3),
"\nP_SNPWISE_MEAN: ", formatC(filebmi$P_SNPWISE_MEAN, format = "e", digits = 3),
"\nP_SNPWISE_TOP1: ", formatC(filebmi$P_SNPWISE_TOP1, format = "e", digits = 3),
"\nN: ", filebmi$N) , NA)
#Run Hudson
igmirror(top=file, bottom=filebmi,
tline=2.75e-06,
bline=2.75e-06,
toptitle="GWAS of Vegetarianism (gene-level; no BMI in model)",
bottomtitle = "GWAS of Vegetarianism (gene-level; BMI adjusted)",
highlight_p = c(1e-3,1e-3),
highlighter="#2DEDDD",
file=paste0(outdir, "/PlotVegGWAS.genelevel")
)pheno1=c("Alanine_aminotransferase_rint", "Albumin_rint", "Alkaline_phosphatase_rint",
"Apolipoprotein_A_rint", "Apolipoprotein_B_rint", "Aspartate_aminotransferase_rint" ,
"C_reactive_protein_rint", "Calcium_rint", "Cholesterol_rint", "Creatinine_rint" ,
"Cystatin_C_rint", "Direct_bilirubin_rint", "Gamma_glutamyltransferase_rint",
"HbA1c_rint", "HDL_cholesterol_rint", "IGF_1_rint" ,
"LDL_direct_rint", "Lipoprotein_A_rint",
"Phosphate_rint", "SHBG_rint", "Testosterone_rint" ,
"Total_bilirubin_rint", "Total_protein_rint", "Triglycerides_rint",
"Urate_rint", "Urea_rint", "Vitamin_D_rint", "eGFR_rint" ,
"bioavailableTest_rint", "freeTest_rint")
phenob=paste0(pheno, ".adjBMI.1df.tsv")
pheno=paste0(pheno1, ".1df.tsv")
phenotypenames=c("Alanine aminotransferase", "Albumin", "Alkaline phosphatase",
"Apolipoprotein A", "Apolipoprotein B", "Aspartate aminotransferase" ,
"C-reactive protein", "Calcium", "Cholesterol", "Creatinine" ,
"Cystatin C", "Direct bilirubin", "Gamma glutamyltransferase",
"HbA1c", "HDL cholesterol", "IGF 1" ,
"LDL direct", "Lipoprotein A",
"Phosphate", "SHBG", "Total testosterone" ,
"Total bilirubin", "Total protein", "Triglycerides",
"Urate", "Urea", "Vitamin D", "eGFR" ,
"Bioavailable testosterone", "Free testosterone")
dir="/Users/mike/Documents/Research/Vegetarian/GEM5/summary/noBMI/"
dirbmi="/Users/mike/Documents/Research/Vegetarian/GEM5/summary/adjBMI/"
outdir="Interactive.GWIS.variant"
for (p in 1:length(pheno)){
file=vroom(paste0(dir, pheno[p]), show_col_types = FALSE,progress = FALSE)
filebmi=vroom(paste0(dirbmi, phenob[p]), show_col_types = FALSE,progress = FALSE)
filebmi
#Set colnames to hudson format
colnames(file)=c("SNP", "CHR", "POS",
"ALLELE0","ALLELE1","N","A1FREQ",
"N_nonveg","effect_allele_frequency_nonveg","N_veg",
"effect_allele_frequency_veg", "beta_marginal" ,
"robust_SE_beta_marginal", "beta_mainG", "beta_interaction",
"robust_SE_beta_mainG","robust_SE_beta_interaction",
"robust_Cov_main_int","robust_p_value_Marginal", "pvalue" ,
"robust_p_value_joint")
colnames(filebmi)=c("SNP", "CHR", "POS",
"ALLELE0","ALLELE1","N","A1FREQ",
"N_nonveg","effect_allele_frequency_nonveg","N_veg",
"effect_allele_frequency_veg", "beta_marginal" ,
"robust_SE_beta_marginal", "beta_mainG", "beta_interaction",
"robust_SE_beta_mainG","robust_SE_beta_interaction",
"robust_Cov_main_int","robust_p_value_Marginal", "pvalue" ,
"robust_p_value_joint")
#Add Link col
file$Link=ifelse(file$pvalue<5e-05,
paste0("https://www.ncbi.nlm.nih.gov/snp/",
file$SNP),
NA)
filebmi$Link=ifelse(filebmi$pvalue<5e-05,
paste0("https://www.ncbi.nlm.nih.gov/snp/",
filebmi$SNP),
NA)
#Add Hover col
file$Hover= ifelse(file$pvalue<5e-05,
paste0("rsID: ", file$SNP,
"\nCHR:POS: ", file$CHR, ":", file$POS,
"\nEA/NEA: ", paste0(file$ALLELE1, "/", file$ALLELE0),
"\nEA freq: ", file$A1FREQ,
"\nBeta (SE) Interaction : ", paste0(signif(file$beta_interaction, digits = 3), " (", signif(file$robust_SE_beta_interaction, digits = 3),")"),
"\nP-value interaction: ", formatC(file$pvalue, format = "e", digits = 3),
"\nBeta (SE) Marginal : ", paste0(signif(file$beta_marginal, digits = 3), " (", signif(file$robust_SE_beta_marginal, digits = 3),")"),
"\nP-value marginal: ", formatC(file$robust_p_value_Marginal, format = "e", digits = 3),
"\nN: ", file$N),
NA)
filebmi$Hover= ifelse(filebmi$pvalue<5e-05,
paste0("rsID: ", filebmi$SNP,
"\nCHR:POS: ", filebmi$CHR, ":", filebmi$POS,
"\nEA/NEA: ", paste0(filebmi$ALLELE1, "/", filebmi$ALLELE0),
"\nEA freq: ", filebmi$A1FREQ,
"\nBeta (SE) Interaction : ", paste0(signif(filebmi$beta_interaction, digits = 3), " (", signif(filebmi$robust_SE_beta_interaction, digits = 3),")"),
"\nP-value interaction: ", formatC(filebmi$pvalue, format = "e", digits = 3),
"\nBeta (SE) Marginal : ", paste0(signif(filebmi$beta_marginal, digits = 3), " (", signif(filebmi$robust_SE_beta_marginal, digits = 3),")"),
"\nP-value marginal: ", formatC(filebmi$robust_p_value_Marginal, format = "e", digits = 3),
"\nN: ", filebmi$N),
NA)
#Run Hudson
igmirror(top=file, bottom=filebmi,
tline=5e-08,
bline=5e-08,
toptitle=paste(phenotypenames[p], "Interaction"),
bottomtitle = paste(phenotypenames[p], "Interaction (BMI-adjusted)"),
highlight_p = c(5e-05, 5e-05),
highlighter="#36E072",
file=paste0(outdir, "/",pheno1[p])
)
}dir="/Users/mike/Documents/Research/Vegetarian/MAGMA5/sumstats/noBMI"
dirbmi="/Users/mike/Documents/Research/Vegetarian/MAGMA5/sumstats/adjBMI"
outdir="Interactive.GWIS.gene"
#Prep gene annotation file
geneloc=as_tibble(read.table(paste0(dir, "/NCBI37.3.gene.loc")))
colnames(geneloc)=c("GENE","CHR","START","STOP", "STRAND", "GENE_NAME")
geneloc=geneloc%>%filter(CHR!="X" & CHR!="Y")
geneloc$CHR=as.integer(geneloc$CHR)
geneloc=geneloc%>%select(GENE, GENE_NAME)
for (p in 1:length(pheno)){
file=read_table(paste0(dir, "/",pheno1[p],"/2.geneanalysis.genes.out"))
filebmi=read_table(paste0(dirbmi, "/",pheno1[p],"/2.geneanalysis.genes.out"))
#only use ID and name as start and stop will be different because of the +2kb and -1kb window added during analysis
file=left_join(file,geneloc)
file=file%>%select(GENE_NAME, everything())
filebmi=left_join(filebmi,geneloc)
filebmi=filebmi%>%select(GENE_NAME, everything())
#Set colnames to hudson format
colnames(file)=c("SNP","GENE_ID_NUMBER", "CHR", "POS",
"STOP","NSNPS","NPARAM",
"N","ZSTAT","pvalue",
"P_SNPWISE_MEAN", "P_SNPWISE_TOP1")
colnames(filebmi)=c("SNP","GENE_ID_NUMBER", "CHR", "POS",
"STOP","NSNPS","NPARAM",
"N","ZSTAT","pvalue",
"P_SNPWISE_MEAN", "P_SNPWISE_TOP1")
#Add Link col
file$Link=ifelse(file$pvalue<1e-3,paste0("https://www.ncbi.nlm.nih.gov/gene/?term=", file$SNP),NA)
filebmi$Link=ifelse(file$pvalue<1e-3,paste0("https://www.ncbi.nlm.nih.gov/gene/?term=", filebmi$SNP),NA)
#Add Hover col
file$Hover= ifelse(file$pvalue<1e-3,paste0("Gene: ", file$SNP,
"\nStart-Stop: ", file$CHR, ":", file$POS, "-", file$STOP,
"\nNSNPS: ", file$NSNPS,
"\nNPARAM: ", file$NPARAM,
"\nP_MULTI: ", formatC(file$pvalue, format = "e", digits = 3),
"\nP_SNPWISE_MEAN: ", formatC(file$P_SNPWISE_MEAN, format = "e", digits = 3),
"\nP_SNPWISE_TOP1: ", formatC(file$P_SNPWISE_TOP1, format = "e", digits = 3),
"\nN: ", file$N), NA)
filebmi$Hover= ifelse(filebmi$pvalue<1e-3,paste0("Gene: ", filebmi$SNP,
"\nStart-Stop: ", filebmi$CHR, ":", filebmi$POS, "-", filebmi$STOP,
"\nNSNPS: ", filebmi$NSNPS,
"\nNPARAM: ", filebmi$NPARAM,
"\nP_MULTI: ", formatC(filebmi$pvalue, format = "e", digits = 3),
"\nP_SNPWISE_MEAN: ", formatC(filebmi$P_SNPWISE_MEAN, format = "e", digits = 3),
"\nP_SNPWISE_TOP1: ", formatC(filebmi$P_SNPWISE_TOP1, format = "e", digits = 3),
"\nN: ", filebmi$N) , NA)
#Run Hudson
igmirror(top=file, bottom=filebmi,
tline=2.75e-06,
bline=2.75e-06,
toptitle=paste(phenotypenames[p], "Interaction gene-level"),
bottomtitle = paste(phenotypenames[p], "Interaction (BMI-adjusted) gene-level"),
highlight_p = c(1e-3,1e-3),
highlighter="#2DEDDD",
file=paste0(outdir, "/", pheno1[p],".gene")
)
}