Medical_code/CAD_autophagy.R

# 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")
# install.packages("tidyverse", repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")

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)


######################################################################
######################################################################



# 本地文件孟德尔分析
# Major coronary heart disease event
# Coronary atherosclerosis
# install.packages("vroom")
#加载包
library(TwoSampleMR)
library(vroom)
library(tidyverse)
library(tools)
library(ieugwasr)
setwd("/root/r_base")

#读取eaf文件
eaf_dat_ori <- vroom("protein_data/eaf_data/assocvariants.annotated.txt.gz", col_names = TRUE)
eaf_dat = eaf_dat_ori[, c(3,7)]
colnames(eaf_dat)

# 指定文件夹路径
folder_path <- "protein_data"

# 获取所有.txt.gz文件
file_list <- list.files(path = folder_path, pattern = "\\.txt\\.gz$", full.names = TRUE)

# 遍历每个文件
for (file in file_list) {
  # 输出当前处理的文件名
  file_name = file_path_sans_ext(file_path_sans_ext(basename(file)))
  #检查是否已经处理过
  if (file.exists(paste0("protein_rds_data/",file_name, ".rds"))) {
    print("文件存在")
    next
  } 
  print(paste("Processing file:", file_path_sans_ext(file_path_sans_ext(basename(file)))))
  #读取暴露
  expo_dat_ori <- vroom(file, col_names = TRUE)
  # 暴露p < 5e-08过滤一下,强相关
  print("=========暴露p < 5e-08过滤一下,强相关=================")
  expo_dat <- subset(expo_dat_ori, expo_dat_ori$Pval<5e-08)

  #eaf填充
  print("=========eaf填充=================")
  expo_dat <- left_join(expo_dat, eaf_dat, by = "Name")
  print(colnames(expo_dat))
  print(dim(expo_dat))

  #标准化列名
  expo_dat = expo_dat[, c(1,2,4,5,6,7,8,10,11,13)]
  colnames(expo_dat) = c("chr", "pos", "SNP", "effect_allele", "other_allele", "beta", "pval", "se", "samplesize", "eaf")
  expo_dat$phenotype_col <- strsplit(file_name, "_")[[1]][3]
  #进行格式转换
  print("=========格式转换=================")
  expo_dat <- TwoSampleMR::format_data(
    expo_dat,
    type = "exposure",
    phenotype_col = "phenotype_col"
  )
  print(colnames(expo_dat))
  print(dim(expo_dat))
  #去除连锁不平衡
  print("=========去除连锁不平衡=================")
  # expo_dat_clump <- clump_data(expo_dat, clump_kb = 10000, clump_r2 = 0.001)
  exposure_data_clumped <- ld_clump(
    # 指定包含`rsid`、`pval`和`id`的数据框
    dat = tibble(rsid = expo_dat$SNP,  # 从 exposure_data 中提取 SNP 标识符
                pval = expo_dat$pval.exposure,  # 从 exposure_data 中提取关联 p 值
                id = expo_dat$exposure),  # 从 exposure_data 中提取 SNP 的 id
    clump_kb = 10000,      # 设置 clumping 窗口大小
    clump_r2 = 0.001,      # 指定用于 LD 的 r^2 阈值
    clump_p = 1,           # 设置 p 值阈值为 1，这样所有的 SNP 都会被纳入 clump，因为咱们前面已经进行过 p 值过滤了嘛
    bfile = "clump_file/EUR",       # 指定 LD 的参考数据集，注意修改你的路径
    plink_bin = "clump_file/plink"  # 使用指定目录下的 plink 二进制文件，注意修改你的路径
  )
  exposure_data_clumped <- expo_dat %>% filter(expo_dat$SNP %in% exposure_data_clumped$rsid)
  print(colnames(exposure_data_clumped))
  print(dim(exposure_data_clumped))
  #保存为rds
  saveRDS(exposure_data_clumped, paste0("protein_rds_data/",file_name, ".rds"))
}

# 读取结局
outcome_dat_ori <- vroom("CHD/finngen_R11_I9_CORATHER.gz", col_names = TRUE)
colnames(outcome_dat_ori)
# 标准化列名
outcome_dat_ori = outcome_dat_ori[, c(1,2,3,4,5,7,9,10,11)]
colnames(outcome_dat_ori) = c("chr", "pos", "other_allele", "effect_allele","SNP", "pval", "beta", "se","eaf")
outcome_dat_ori$phenotype_col <- "finngen_R11_I9_CORATHER"

# 循环进行等位基因对齐
# 指定文件夹路径
folder_rds_path <- "protein_rds_data"

# 获取所有.txt.gz文件
file_rds_list <- list.files(path = folder_rds_path, pattern = "\\.rds$", full.names = TRUE)
# 遍历每个文件进行等位基因对齐
for (file_rds in file_rds_list) {
  # 输出当前处理的文件名
  file_name = file_path_sans_ext(file_path_sans_ext(basename(file_rds)))
  print(paste("Processing file:", file_path_sans_ext(file_path_sans_ext(basename(file_rds)))))
  #读取rds
  exposure_data_clumped = readRDS(file_rds)
  #暴露和结局取交集，确定哪些SNP是两者共有的
  print("=========暴露和结局取交集=================")
  data_common <- merge(exposure_data_clumped, outcome_dat_ori, by = "SNP")
  print(dim(data_common))
  print(colnames(data_common))
  if (nrow(data_common) > 0) {
    # 提取结局数据，并进行格式修改，以用于后续的 MR 分析
    outcome_dat <- format_data(outcome_dat_ori,
                                # 设置类型为 "outcome"，表示这是一个结局数据
                                type = "outcome",
                                # 选择用于 MR 分析的 SNP，这里使用了 data_common 数据集中的 SNP 列
                                snps = data_common$SNP,
                                # 设置 SNP 列的名称
                                snp_col = "SNP",
                                # 设置效应大小 (beta) 的列名称
                                beta_col = "beta",
                                # 设置标准误差 (SE) 的列名称
                                se_col = "se",
                                # 设置效应等位基因频率 (EAF) 的列名称
                                eaf_col = "eaf",
                                # 设置效应等位基因的列名称
                                effect_allele_col = "effect_allele",
                                # 设置非效应等位基因的列名称
                                other_allele_col = "other_allele",
                                # 设置 p 值 (p-value) 的列名称
                                pval_col = "pval",
                                # 设置样本大小 (sample size) 的列名称
                                samplesize_col = "samplesize",
                                # 设置病例数 (number of cases) 的列名称
                                ncase_col = "ncase",
                                # 设置 id 的列名称
                                id_col = "id",
                                phenotype_col = "phenotype_col")

    #等位基因对齐
    print("=========等位基因对齐=================")
    dat <- harmonise_data(exposure_dat = exposure_data_clumped, outcome_dat = outcome_dat)
    print(dim(dat))
    print(colnames(dat))
    #保存为rds
    saveRDS(dat, paste0("harmonise_data_CORATHER/",file_name, ".rds"))
  } else {
    print("数据条数为 0")
  }
}

# 循环进行MR分析
# 指定文件夹路径
folder_rds_path_m <- "harmonise_data_CHD"

# 获取所有.txt.gz文件
file_rds_list_m <- list.files(path = folder_rds_path_m, pattern = "\\.rds$", full.names = TRUE)

# 遍历每个文件进行mr分析
result_all <- data.frame()
for (file_rds_m in file_rds_list_m) {
  #读取rds
  dat = readRDS(file_rds_m)
  #去掉结局相关SNP
  dat <- subset(dat, dat$pval.outcome>=5e-08)
  print("=========MR 分析=================")
  # MR 分析
  result <- mr(dat)
  result_all <- rbind(result_all, result)
}
# 保存为 CSV 文件
write.csv(result_all, file = "result_all_CHD.csv", row.names = FALSE)

a = read.csv("result_all_CHD.csv")
View(a)

# b = read.csv("result_all_CHD.csv")
# View(b)


#单独跑MR分析
dat = readRDS("harmonise_data_CORATHER_result/18387_7_ST13_suppression_of_tumorigenicity_13.rds")
dat <- subset(dat, dat$pval.outcome>=5e-08)
result <- mr(dat)
View(result)
#生成OR值 
generate_odds_ratios(result)

# 异质性检验
mr_heterogeneity(dat)

# 水平多效性检验
mr_pleiotropy_test(dat)

# 散点图
p1 <- mr_scatter_plot(result, dat)
p1

# 森林图
result_single <- mr_singlesnp(dat)
p2 <- mr_forest_plot(result_single)
p2

# 留一图
result_loo <- mr_leaveoneout(dat)
p3 <- mr_leaveoneout_plot(result_loo)
p3

# 漏斗图
result_single <- mr_singlesnp(dat)
p4 <- mr_funnel_plot(result_single)
p4