解螺旋孟德尔随机化实操代码练习记录

## 解螺旋孟德尔随机化9-1代码练习记录

##一、准备工作，安装R包

install.packages('devtools')
devtools::install_github('MRCIEU/TwoSampleMR')

library(TwoSampleMR) #加载R包

#获取谷歌授权的IEU OPEN GWAS API的vpn：https://api.opengwas.io/
Sys.setenv(OPENGWAS_JWT="TOKEN编码")



##二、导入暴露数据
#利用IEU数据库导入暴露数据
#下述函数代码已包括了P值的筛选及clump连锁不平衡的筛选，默认是P<5x10-8，若需更改p值，则在clump前加“p1=1e-5”;clump函数是对LD的SNP集进行修剪，否则使用相关的SNP会导致重复计算，从而导致有偏倚的因果效应估计。因此通过设置r2和kb参数来制定去除LD的标准，默认的是r2=0.001，kb=10000。accesstoken这个参数，对于中国大陆地区的用户必须设置该参数为accesstoken=NULL，这样才能顺利获取数据，否则就需要开VPN获取谷歌授权。
CHE_exposure<-extract_instruments (outcomes= 'ukb-b-1489',clump=TRUE)



##三、导入结局(获取暴露因素-结局变量的GWAS数据)
#SNP选择的是暴露因素中选取的SNPS，也就是上面的65个，outcomes就是在IEU数据库中的ID，可以直接复制过来。
#proxies表示是否使用代理SNP，默认值是TRUE，也即当一个SNP在outcome中找不到时可以使用与其存在强连锁不平衡的SNP信息来替代，看个人喜好设置为FALSE或TURE。
#maf threshold表示的是SNP在outcome中的最小等位基因频率，默认值是0.3，不过大样本GWAS可以适当调低，这里设置的是0.01。
#access token的意思同上，大陆用户必须设置成access token=NULL。
CHE_data<-extract_outcome_data(
  snps = CHE_exposure$SNP,
  outcomes ="ebi-a-GCST009541",
  proxies = FALSE,
  maf_threshold =0.01)



##四、对数据进行预处理，使其效应等位与效应量保持统一
#参数action最重要，一般我们推荐使用默认值action=2即可，当然也可以使用action=3，这时候就表示去除所有存在回文结构的SNP。
CHE_H_data<-harmonise_data(
  exposure_dat=CHE_exposure,
  outcome_dat =CHE_data,
  action=2)



##五、MR分析
#b值代表效应值的意思，若是负值，则代表暴露因素的增加会导致结局变量减少，即奶酪摄入增加会导致心衰的风险降低，和原文的结论一致。
#此外，如果针对的是二分类结局那么需要把B转化为风险比OR，而OR值则是大于1增加风险，小于1则降低风险.
#P值代表我们的暴露因素与结局变量是否具有统计学意义。
mr_CHE_results<-mr(CHE_H_data)
mr_CHE_results
OR_CHE_results<- generate_odds_ratios(mr_CHE_results)
OR_CHE_results



##六、敏感性分析
het_CHE<-mr_heterogeneity(CHE_H_data)#异质性检验
#该代码默认的是固定效应模型。但通过分析可以看到，这些IV之间存在很强的异质性(Qpva1远小于0.05)，这时候我们需要剔除某些outcome的P值非常小的SNP.
#或者直接使用随机效应模型来估计MR效应量(即MR-ivw分析的结果)
het2_CHE<- mr(CHE_H_data,method_list=c('mr_ivw_mre'))

pleio_CHE<-mr_pleiotropy_test(CHE_H_data)#水平多效性检验
pleio_CHE



##七、结果可视化#散点图
CHE_plotl<-mr_scatter_plot(mr_CHE_results, CHE_H_data)
CHE_plotl

#逐个剔除检验，留一法
CHE_res_loo<-mr_leaveoneout(CHE_H_data)
CHE_plot2<-mr_leaveoneout_plot(CHE_res_loo)
CHE_plot2

#单个snp获得的估计来生成森林图
CHE_res_single<-mr_singlesnp(CHE_H_data)
CHE_plot3<-mr_forest_plot(CHE_res_single)
CHE_plot3

#漏斗图
CHE_plot4<-mr_funnel_plot(CHE_res_single)
CHE_plot4






## 解螺旋孟德尔随机化9-2代码练习记录


##一、准备工作，安装R包

install.packages('devtools')
devtools::install_github('MRCIEU/TwoSampleMR')

library(TwoSampleMR) #加载R包

#获取谷歌授权的IEU OPEN GWAS API的vpn：https://api.opengwas.io/
Sys.setenv(OPENGWAS_JWT="TOKEN编码")


#获取暴露数据
che<-extract_instruments(outcomes='ukb-b-1489', clump=TRUE)#获取奶酪摄入snp
che


#剔除混杂因素 http://www.phenoscanner.medsch1.cam.ac.uk/   #网站好卡，还没成功进去过，只能换时间段继续尝试
#将获取的snp号（如：rs78876700），一个一个放进网站搜索，看该snp的表型是否包含我们查找出来的混杂因素

#视频发现第一行snp包含混杂因素，使用以下函数去除第1行snp
che<-che[-1,]#删除某行
che
write.csv(che,file='che.csv')#把整理后的snp表保存成csv文件
write.csv
View(write.csv)


#获取相应的结局数据 #冠心病数据
#snps：它是一串以rs开头的SNP ID
#outcomes：它是outcome在MR base中的ID；
#proxies：它表示是否使用代理SNP，默认值是TRUE，也即当一个SNP在outcome中找不到时可以使用与其存在强连锁不平衡的SNP信息来替代，我个人喜欢设置成FALSE。
#maf_threshold：它表示的是SNP在outcome中的最小等位基因频率，默认值是0.3，不过大样本GWAS可以适当调低，这里设置的是0.01。
chd<-extract_outcome_data(
  snps = che$SNP,
  outcomes ='ieu-a-7',
  proxies = FALSE,
  maf_threshold = 0.01)
#找出的结局数据的snp数量一定少于或等于暴露数据，因为这一步是用暴露数据到结局数据中进行匹配、提取


#把暴露数据和结局数据的效应等位基因对齐（将IV的效应等位基因（effect allele）对齐），对齐后才能进行MR分析
#这一步代码一般不用改
#对齐过程中，可能会去掉一些无法对齐的数据，注意留意数据数量有无改变
mydata <- harmonise_data(
  exposure_dat=che,
  outcome_dat=chd,
  action= 2)
mydata


#MR分析
res <- mr(mydata)
res
OR<-generate_odds_ratios(res)
OR

#异质性检验
het <- mr_heterogeneity(mydata)
het
#结果：这些IV之间存在很强的异质性（Q_pval远小于0.05），
#这时候我们需要剔除某些outcome的P值非常小的SNP，
#或者直接使用随机效应模型来估计MR效应量
mr(mydata,method_list=c('mr_ivw_mre'))#随机效应模型

#多效性检验
pleio <- mr_pleiotropy_test(mydata)
pleio

#逐个剔除检验
single <- mr_leaveoneout(mydata)
mr_leaveoneout_plot(single)

#散点图
mr_scatter_plot(res,mydata)

#森林图
res_single <- mr_singlesnp(mydata)
mr_forest_plot(res_single)

#漏斗图
mr_funnel_plot(res_single)