# apolipoprotein L, 1
# apolipoprotein L, 1
# arylsulfatase B
# Baculoviral IAP repeat-containing protein 5
# BCL2/adenovirus E1B 19 kDa protein-interacting protein 3
# calnexin
# Caspase 8 levels
# Caspase-3
# C-C motif chemokine ligand 2
# Cyclin-dependent kinase inhibitor 1B
# Cathepsin L1 levels
# cathepsin D
# Cathepsin B
# Cathepsin L1 levels
# C-X3-C motif chemokine ligand 1
# Death-associated protein kinase 2
# DnaJ homolog subfamily B member 9
# Glutamate receptor ionotropic, delta-2
# Interferon gamma
# Gamma-aminobutyric acid receptor-associated protein-like 1
# Gamma-aminobutyric acid receptor-associated protein-like 2
# Glutamate receptor ionotropic, delta-2
# Mitogen-activated protein kinase 1
# Mitogen-activated protein kinase 3
# Mitogen-activated protein kinase 8
# Nicotinamide phosphoribosyltransferase
# Next to BRCA1 gene 1 protein
# PTK6
# sphingosine kinase 1
# Vesicle-associated membrane protein 3
# vascular endothelial growth factor A
#在线孟德尔分析
#加载包
library(TwoSampleMR)
# install.packages("ieugwasr")
library(ieugwasr)
# 登陆token
Sys.setenv(OPENGWAS_JWT="eyJhbGciOiJSUzI1NiIsImtpZCI6ImFwaS1qd3QiLCJ0eXAiOiJKV1QifQ.eyJpc3MiOiJhcGkub3Blbmd3YXMuaW8iLCJhdWQiOiJhcGkub3Blbmd3YXMuaW8iLCJzdWIiOiJ6anQ3ODU2MzIxQGhvdG1haWwuY29tIiwiaWF0IjoxNzI0NTg2NzM3LCJleHAiOjE3MjU3OTYzMzd9.hseY1Y4YiYJ4WWi-WoQmm0VWQA6Ozqurj08HvJV8q6xh7D4Jfv5o3zdwBvU_4HhC87C59KWgXaBm5ZiMynqLSp5uFLlai_5ObswpbdLHtNRghZYti0wF7TRTGcXPa9_ICWvrEoIRNoAyj10I1nOrPE9RhRL0PrzHOYb9Y3L7oBDzmsch-boFLYE3R9ZhUoZSMMiZYNR9_PmYbJi2Lrd6C_tI1fCna9jzF5vHsQQReYgy4IS59rr_lDEwdE46dCIq2cHoDIZJ4ohArPHtoLaw5CXf0yUrY89Xo2J5KnVWHtzlFK7EEgOIpx1DX6mSgrVxXtvlFSM-j7K1qbM9Awr6JQ")

ieugwasr::user()

Coronary_artery_disease_list = c("ebi-a-GCST005195", "ebi-a-GCST005194")

autoghagy_list <- c("prot-a-134", "prot-a-135", "prot-a-172", "prot-a-252" , "prot-a-261",
            "prot-a-353", "ebi-a-GCST90012072", "prot-a-364", "prot-b-4", "prot-a-495",
            "ebi-a-GCST90012073", "prot-b-51", "prot-a-718", "ebi-a-GCST90012073", "prot-b-70",
            "prot-a-763", "prot-a-842", "prot-a-1276", "prot-a-1428", "prot-a-1161", 
            "prot-a-1162", "prot-a-1276","prot-a-1845", "prot-a-1848", "prot-a-1849", 
            "prot-a-1998", "prot-a-2004","prot-a-2536","ebi-a-GCST90000502")

result_all <- data.frame()
for (i in seq_along(Coronary_artery_disease_list)) {
  for (j in seq_along(autoghagy_list)) {
    #两样本MR(在线读)
    #提取蛋白的SNP
    print(autoghagy_list[j])
    autoghagy <- extract_instruments(outcomes = autoghagy_list[j])
    #提取SNP在结局中的信息
    print(Coronary_artery_disease_list[i])
    outcome_dat <- extract_outcome_data(snps = autoghagy$SNP, outcomes = Coronary_artery_disease_list[i])
    #合并暴露和结局的数据
    dat <- harmonise_data(autoghagy, outcome_dat)
    #分析MR
    result <- mr(dat)
    result_all <- rbind(result_all, result)
  }
}
# 保存为 CSV 文件
write.csv(result_all, file = "result_all.csv", row.names = FALSE)

a = read.csv("result_all.csv")

autoghagy <- extract_instruments(outcomes = "ebi-a-GCST90000502")
#提取SNP在结局中的信息
print(Coronary_artery_disease_list[i])
outcome_dat <- extract_outcome_data(snps = autoghagy$SNP, outcomes = "ebi-a-GCST005194")
#合并暴露和结局的数据
dat <- harmonise_data(autoghagy, outcome_dat)
#分析MR
result <- mr(dat)

# #由于结局是二分类，所以生成OR值
# OR <- generate_odds_ratios(results)

# # 异质性检验
# heterogeneity <- mr_heterogeneity(dat)
# res_MRPRESSO <- run_mr_presso(dat)

# # 多效性检验
# pleiotropy <- mr_pleiotropy_test(dat)
# pleiotropy

# #散点图
# mr_scatter_plot(results, dat)

# #留一图
# leaveoneout <- mr_leaveoneout(dat)
# mr_leaveoneout_plot(leaveoneout)

# #森林图
# results_single <- mr_singlesnp(dat)
# mr_forest_plot(results_single)

# #漏斗图
# mr_funnel_plot(results_single)

# 本地文件孟德尔分析
install.packages("vroom")
library(vroom)
setwd("/root/zjt")
a <- vroom("11510_31_APOL1_Apo_L1.txt", col_names = TRUE)