rm(list=ls())
getwd()
setwd( "E:/Rproject/05_Dental caries_AD_MR/Periodontal disease/PEERJ")

#----1. packages installation----
library(TwoSampleMR) 
library(MRPRESSO)
library(forestplot)
library(Ipaper)



#-------DENTAL CARIES ON AD------
#----2. exposure profiles----
## 1.1.exposure data----
exp.local <- read_exposure_data(
  filename = 'dental_caries_47IVs_GWAS.csv',
  clump = FALSE,
  samplesize = 'N',
  sep= ",",
  snp_col = "SNP",
  beta_col = "BETA",
  se_col = "SE",
  effect_allele_col ="EA",
  other_allele_col = "OA",
  eaf_col = "EAF",
  pval_col = "Pvalue",
  phenotype_col = "traits",
  id_col = "id",
)
exp <- exp.local

## 1.2 clump_data----
exp <- clump_data(exp,clump_kb=10000,clump_r2=0.001)

## 1.3 remove confounders----
confounders <- c('rs3865314','rs11676272','rs62106258',
                 'rs9366651','rs898797','rs10851907',
                 'rs8054556','rs28822480','rs7429279') #BMI and smoking
exp <- exp[!exp$SNP %in% confounders,]

save(exp,file='exp_DentalCaries.Rdata')
load('exp_DentalCaries.Rdata')

#----3. outcome profiles----
## 3.1 outcome dataframe----
gwasID_outcomes <- c('ieu-b-2', #AD 2019
                     'finn-b-AD_EO',  #Alzheimer disease (Early onset) 
                     'finn-b-AD_LO',  #Alzheimer disease (Late onset)
                     'finn-b-G6_ALZHEIMER' #Alzheimer disease FInn
)
gwasID_m.out <- matrix(gwasID_outcomes,1,4) 
rownames(gwasID_m.out) <- 'gwasID_outcomes'
colnames(gwasID_m.out) <- c('AD19','ADEO','ADLO','FinnAD')
gwasID_dat.out <- data.frame(gwasID_m.out)

resnames <- matrix(nrow=1,ncol=ncol(gwasID_dat.out))
for(i in 1:(ncol(gwasID_dat.out))){
  resnames[1,i] <- paste('DentalCaries',colnames(gwasID_dat.out)[i],sep="_")
}

## 3.2 extract outcome data----
out_subtype <- list()
out <- paste('DentalCaries',colnames(gwasID_dat.out),sep='_')
for(j in 1:(ncol(gwasID_dat.out))){
  out_subtype[[j]] = extract_outcome_data(snps=exp$SNP,
                                          outcomes= gwasID_dat.out [1,j],
                                          proxies=TRUE,
                                          maf_threshold=0.01,
                                          access_token=NULL)
}
save(gwasID_dat.out,out_subtype,file='out_AD.Rdata')
load('out_AD.Rdata')

#----4. harmonise data----
data=list()
for(j in 1:(ncol(gwasID_dat.out))){
  data[[j]] <- harmonise_data(exp,out_subtype[[j]],action=2)
}

#----5. MR analysis----
res=list()
for(j in 1:(ncol(gwasID_dat.out))){res[[j]] <- mr(data[[j]],
                                                  method_list=c('mr_ivw',
                                                                'mr_egger_regression',
                                                                'mr_simple_median',
                                                                'mr_weighted_median'
                                                  ))}
resOR <- list()
for(j in 1:(ncol(gwasID_dat.out))){resOR[[j]]<- generate_odds_ratios(res[[j]])} #得到相应系数的OR和置信区间
resOR[[1]] #AD19
resOR[[2]] #ADEO
resOR[[3]] #ADLO
resOR[[4]] #FinAD


#----6. sensitivity analysis----
## 6.1 heterogeneity----
het=list()
for(j in 1:(ncol(gwasID_dat.out))){
  het[[j]] <- mr_heterogeneity(data[[j]])
}
## 6.2 pleiotropy-----
### 1). egger
pleio=list()
for(j in 1:(ncol(gwasID_dat.out))){
  pleio[[j]] <- mr_pleiotropy_test(data[[j]])
}
### 2) presso
presso=list()
for(j in 1:(ncol(gwasID_dat.out))){
  presso[[j]] <-mr_presso(BetaOutcome = "beta.outcome",
                          BetaExposure = "beta.exposure",
                          SdOutcome = "se.outcome",
                          SdExposure = "se.exposure",
                          OUTLIERtest = TRUE,
                          DISTORTIONtest = TRUE,
                          data = data[[j]],
                          NbDistribution = 2000,
                          SignifThreshold = 0.05)
}

## 6.3 leave_one_out-----
single=list()
for(j in 1:(ncol(gwasID_dat.out))){
  single[[j]] <- mr_leaveoneout(data[[j]])
}
leaveoneout=list()
for(j in 1:(ncol(gwasID_dat.out))){
  leaveoneout[[j]] <- mr_leaveoneout_plot(single[[j]])
}

save(resOR,het,pleio,presso,single,leaveoneout,
     file='DentalCaries_AD_MRresults.Rdata')

#----7. visulization----
## 7.1 scatter plots----
p1_scatter=list()
for(j in 1:(ncol(gwasID_dat.out))){
  p1_scatter[[j]] <- mr_scatter_plot(res[[j]],data[[j]])
}
## 7.2 forest plots  single----
res_single=list()
for(j in 1:(ncol(gwasID_dat.out))){
  res_single[[j]] <- mr_singlesnp(data[[j]])
}
p1_forest=list()
for(j in 1:(ncol(gwasID_dat.out))){
  p1_forest[[j]] <- mr_forest_plot(res_single[[j]])
}
## 7.3 leave_one_out plot----
single=list()
for(j in 1:(ncol(gwasID_dat.out))){
  single[[j]] <- mr_leaveoneout(data[[j]])
}
p1_leaveoneout=list()
for(j in 1:(ncol(gwasID_dat.out))){
  p1_leaveoneout[[j]] <- mr_leaveoneout_plot(single[[j]])
}

## 7.4 funnel plot----
p1_funnel=list()
for(j in 1:(ncol(gwasID_dat.out))){
  p1_funnel[[j]] <- mr_funnel_plot(res_single[[j]])
}

## 7.5 save----
save(p1_scatter,p1_forest,p1_leaveoneout,p1_funnel,file='DenralCaries_AD_MR_Figures.Rdata')


#-------AD on DENTAL CARIES------
traits <- c('AD19','FinnAD')
ID <- c('ieu-b-2','finn-b-G6_ALZHEIMER')
traits.df <- as.data.frame(traits,ID)
traits.df$ID <- rownames(traits.df)
#----2. exposure profiles----
exp0 <- list()
for(j in 1:2){
  exp0[[j]] = extract_instruments(outcomes=traits.df[j,2],
                                  p1 = 5e-8, 
                                  clump = TRUE, 
                                  r2 = 0.001,
                                  kb = 10000,
                                  access_token = NULL) 
}
exp <- exp0

#1.3 confounders
##1.3.1 phenoscanner
library(readr)
exp_AD_phenoscanner <- read_tsv('exp_AD_phenoscanner_GWAS.tsv.gz')
exp_AD_phenoscanner


save(exp,file='exp_AD.Rdata')

#----3. outcome profiles----
## 3.1 extract outcome data----
out <- list()    
for(j in (1:2)){
  out[[j]] = extract_outcome_data(snps=exp[[j]]$SNP,
                                  outcomes= 'ukb-b-4770',  #Dental caries, unspecified
                                  proxies=TRUE,
                                  maf_threshold=0.01,
                                  access_token=NULL)
}

save(out,file='out_DentalCaries.Rdata')


#----4. harmonise data----
data=list()
for(j in 1:2){
  data[[j]] <- harmonise_data(exp[[j]],out[[j]],action=2)
}

#----5. MR analysis----
res=list()
for(j in 1:2){res[[j]] <- mr(data[[j]],
                             method_list=c('mr_ivw',
                                           'mr_egger_regression',
                                           'mr_simple_median',
                                           'mr_weighted_median'
                             ))}
res[[1]]
res[[2]]

#----6. sensitivity analysis----
## 6.1 heterogeneity----
het=list()
for(j in 1:2){
  het[[j]] <- mr_heterogeneity(data[[j]])
}


## 6.2 pleiotropy----
### 1) egger----
pleio=list()
for(j in 1:2){
  pleio[[j]] <- mr_pleiotropy_test(data[[j]])
}
### 2) presso----
presso=list()
for(j in 1:2){
  presso[[j]] <-mr_presso(BetaOutcome = "beta.outcome",
                          BetaExposure = "beta.exposure",
                          SdOutcome = "se.outcome",
                          SdExposure = "se.exposure",
                          OUTLIERtest = TRUE,
                          DISTORTIONtest = TRUE,
                          data = data[[j]],
                          NbDistribution = 2000,
                          SignifThreshold = 0.05)
}

## 6.3 leave_one_out----
single=list()
for(j in 1:2){
  single[[j]] <- mr_leaveoneout(data[[j]])
}
leaveoneout=list()
for(j in 1:2){
  leaveoneout[[j]] <- mr_leaveoneout_plot(single[[j]])
}


save(resOR,het,pleio,single,leaveoneout,file='AD_DentalCaries_MRresults.Rdata')

#----7. visulization----
## 7.1 scatter plots----
p_scatter=list()
for(j in 1:2){
  p_scatter[[j]] <- mr_scatter_plot(res[[j]],data[[j]])
}
## 7.2 forest plots  single----
res_single=list()
for(j in 1:2){
  res_single[[j]] <- mr_singlesnp(data[[j]])
}
p_forest=list()
for(j in 1:2){
  p_forest[[j]] <- mr_forest_plot(res_single[[j]])
}
## 7.3 leave_one_out plot----
single=list()
for(j in 1:2){
  single[[j]] <- mr_leaveoneout(data[[j]])
}
p_leaveoneout=list()
for(j in 1:2){
  p_leaveoneout[[j]] <- mr_leaveoneout_plot(single[[j]])
}

## 7.4 funnel plot----
p_funnel=list()
for(j in 1:2){
  p_funnel[[j]] <- mr_funnel_plot(res_single[[j]])
}

save(p_scatter,p_forest,p_leaveoneout,p_funnel,
     file='AD_DentalCaries_MR_Figures.Rdata')