Druggable_genome-wide_MR代码

#分析前准备

# 使用“1s()”函数获取当前环境中的所有变量、函数和对象的名称，然后传递给”list”参数

#将“list”参数传递给“rm()”函数，表示制除这些变量、函数和对象

# 因此，“rm(list=ls())”表示删除当前环境中的所有变量、函数和对象，相当于清空当前环境

rm(list = ls())

#判断是否已经安装”pacman”，如果没有则安装它

if(!require(“pacman”)) install.packages(“pacman”, update=F, ask=F)

#设置Bioconductor镜像地址为中国科技大学的镜像

options(BioC_mirror=”https://mirrors.ustc.edu.cn/bioc/”)

#加载”pacman”包，方便加载其他得到R包

library(“pacman”)

#使用”P_load”函数加载所需要的R包,

#p_load会判断环境中是否有这个包，如果没有，先安装再加载；如果有，那么直接加载。

#第一次需要安装，以后就无须安装，把前面的#去掉

#install.packages(“devtools”)

#unlink(“D:/R-4.2.3”, recursive = TRUE) #Permission denied的问题

#devtools::install_github(“MRCIEU/TwoSampleMR”,force = TRUE)

p_load(ggplot2,readr)

p_load_gh(“MRCIEU/TwoSampleMR”)

#mrcieu/gwasglue”.

#install.packages(“Matrix”)    #如果缺失，则安装此R包

#install.packages(“MendelianRandomization”)  #如果缺失，则安装此R包

library (TwoSampleMR)       #引用包

#暴露数据

#Saturated fatty acids ll id:met-d-SFA

#结局数据 #疾病

#heart failure (HF) id: ebi-a-GCST009541

#暴露数据命名为expofile

expofile=”met-d-Total_FA”

#靶基因位点

#Acat2 160183077..160200144

chr_pos = 6             #染色体位置

pos_start = 160183077   #开始位置

pos_end = 160200144      #结束位置

#结局数据命名为outcfile

outcfile= “ebi-a-GCST009541”

#Step.1 暴露数据处理，在线读取gwas数据

#提取检测的SNP,存到biof_exposure_dat中

biof_exposure_dat <- extract_instruments(outcome= ‘met-d-SFA’)   #”里的为暴露数据

#进行连锁不平衡剔除

biof_exposure_dat <- clump_data(biof_exposure_dat,

              clump_kb=100,      # 定义连锁不平衡窗口大小

              clump_r2=0.3,      # 定义连锁不平衡的R平方阔值

              clump_p1= 5e-08,   # 保留p值最小的SNP

              clump_p2=5e-08)    #剔除p值大于阈值的SNP

#获取去除连锁不平衡后的行数和列数

dim(biof_exposure_dat)

View(biof_exposure_dat)

#Step.2 提取药物靶点周围的SNP

# subset函数从biof exposure_dat筛选SNP

# chr.exposure == chr_pos: chromosome与药物靶点chromosome相等

# pos.exposure >pos_start – 1000000: SNP位置大于把点start位置左侧1000kb

# pos.exposure< pos_end + 1000000: SNP位置小于把点end位置右侧1000kb

#这样提取出药物靶点周围±100kb的SNP

Drug_Target_SNP <- subset(biof_exposure_dat,

              chr.exposure == chr_pos &

                pos.exposure > 160183077 – 1000000 &

                pos.exposure < 160200144 + 1000000)

#将结果写入csv文件

#Target_SNP:提取的药物靶点周围SNP

#“Drug_Target_SNP,csv”:输出的csv文件名

write.csv(Drug_Target_SNP, file =”Drug_Target_SNP.csv”)

View(Drug_Target_SNP)

#Step.3 获取结局中与药物相关SNP

#heart failure (HF) id: ebi-a-GCST90018806

#load outcome data 找到和结局相关性的snp

# 从Outcome数据中提取与药物点SNP相关的表型数据

biof_Outcome_dat <- extract_outcome_data(

  snps= Drug_Target_SNP$SNP,  # 药物靶点SNP

  outcomes = outcfile)        #表型数据文件

# harmonize and merge 数据

#确保SNP对暴露和结果的效应基于同一等位基因

harmonise_dat <- harmonise_data(

  exposure_dat = Drug_Target_SNP,     # 药物靶点SNP数据

  outcome_dat = biof_Outcome_dat,     # 相关表型数据

  action = 1)                         # 调和方法

#查看harmonise后的数据

View(harmonise_dat)

#将结果写入文件

write.table(harmonise_dat,

    “expo_outc_harmonise.csv”,

    row.names = F,

    sep= “\t”,

    quote= F)

#去除混杂因素

p_load(MendelianRandomization.purrr,readr)

#设置变量grp_size,用于设置每个分组的SNP数量

grp_size<- 100

#将data中的SNP分成多个大小为grp_size的分组，保存到变量grps中

grps <- split (harmonise_dat$SNP,ceiling(seq_along(harmonise_dat$SNP)/grp_size))

#通过phenoscanner API对每个分组进行关联分析

#map_dfr将结果合并为一个数据框,map为purrr中的函数，需要先打开

library(purrr)

#对grps中的每个分组运行phenoscanner API进行关联分析

results = map_dfr(grps, ~phenoscanner(snpquery = .x,         # .X表示传入的分组

      catalogue =”GWAS”,       # 在GWAS目录中提索

      pvalue = 1e-05,           # p值闹值设置为le-05

      proxies = “None”,         # 不使用代理SNP

      r2 = 0.8,                 # 相关性阔值设置为0.8

      build = 37) $results)     # 基因组版本为b37

# 将关联结果写入文件confounder07.csv

write_csv(results,”confounder.csv”)

View(results)

# 设置变量remove_snps,用于保存要移除的混杂因素snp

#r如果有的话，更换下方的rs，没有不需要运行

remove_snps<- c(“rs117733303”)

#移除含有混杂SNP的行

harmonise_dat < harmonise_dat[!harmonise_dat$SNP %in% remove_snps,]

# 将移除混杂因素后的结果写入clean_confounder07.csv

write_csv(harmonise_dat,”clean_confounder.csv”)

res=mr(harmonise_dat)

res

write.csv(res,file=”MRres.csv”)

#mendelian randomization (HR) analysis

MR_result<- mr(harmonise_dat)

View(MR_result)

result_OR=generate_odds_ratios(MR_result)

result_OR$or1=1/result_OR$or  #(如果是抑制剂(看药物对靶点的作用)，就用这个倒数，如果不是则删去1/)

#OR值小于1，说明代表暴露因素（脂肪酸）是结局（心衰）的有利因素

#OR值大于1，说明代表暴露因素（脂肪酸）是结局（心衰）的不利因素

result_OR$or_lci951=result_OR$or1-1.96*result_OR$se

result_OR$or_uci951=result_OR$or1+1.96*result_OR$se

write.table(result_OR[,5:ncol(result_OR)],”MR OR.xls”,row.names = F,sep =”\t”,quote = F)

#异质性检验,通过IVW和MR-Egger检验评估异质性，pvalue<0.05说明研究中存在异质性

mr_heterogeneity (harmonise_dat,method_list=c(“mr_egger_regression”, “mr_ivw”))  #异质性检验

#outlier test

run_mr_presso(harmonise_dat,NbDistribution = 1000)  #偏倚检验,snp少的做不出来，没关系

#多效性检验–MR egger，看p值,如果pvalue<0.05，说明数据存在多效性。

#如果存在多效性就需要重新选择工具变量或者重新选择暴露和结局

pleio <- mr_pleiotropy_test(harmonise_dat) 

View (pleio)

write.csv(pleio, file=”pleio.csv”)   #多效性结果

single <- mr_leaveoneout (harmonise_dat)

mr_leaveoneout_plot(single) #留一法检验敏感性

write.csv(single,file=”single.csv”)

p1 <-mr_scatter_plot(res,harmonise_dat)#散点图

p1

ggsave(p1[[1]], file=”scatter.pdf”, width=8, height=8)

#单SNP分析

res_single=mr_singlesnp(harmonise_dat)

#绘制森林图，森林图可以表示每个SNP与结果的关联效

p2 <- mr_forest_plot(res_single)

p2[[1]]

ggsave(p2[[1]], file=”forest.pdf”, width=8, height=8)

#绘制漏斗图

p3 <-mr_funnel_plot(singlesnp_results = res_single)

p3[[1]]

ggsave(p3[[1]], file=”funnel.pdf”, width=8, height=8)

#留一法敏感性分析

p4 <- mr_leaveoneout_plot(leaveoneout_results = mr_leaveoneout(harmonise_dat))

p4[[1]]

ggsave(p4[[1]], file=”leaveoneout.pdf”, width=8, height=8)