学习笔记-孟德尔随机化代码

孟德尔随机化-----本地数据

#打开R包
library(TwoSampleMR)

####1、整理暴露GWAS数据####
setwd("D:/Software/R/projects/MR")#设置工作路径
#读取exposure数据 
a<-read.table("SNP_gwas_mc_merge_nogc.tbl.txt",header = T)
View(a) 

#相关性设置，并将文件放到TwoSampleMR包所在文件夹
b<-subset(a,p<5e-08) 
write.csv(b, file="exposure.csv") 

#将exposure数据读入TwoSampleMR
bmi<-system.file("1- MRlearning/exposure.csv",package = "TwoSampleMR")
bmi_exp_dat<-read_exposure_data(filename = bmi,
								sep = ",",
                                snp_col = "SNP",
                                beta_col = "beta",
                                se_col = "se",
                                effect_allele_col = "effect_allele",
                                other_allele_col = "other_allele",
                                eaf_col = "eaf",
                                clump = FALSE)

#独立性设置，进行clump去除连锁
bmi_exp_dat_clumped<-clump_data(  bmi_exp_dat,
                                  clump_kb = 10000,
                                  clump_r2 = 0.001,
                                  clump_p1 = 1,
                                  clump_p2 = 1,
                                  pop = "EUR")

####2、整理结局GWAS数据####
#读取Alzheimer数据
setwd("D:/Alzheimer IGAP/IGAP_summary_statistics")#设置工作路径
c<-read.table("IGAP_stage_1.txt",header = T)

#将暴露SNP从结局中提取出来，取交集
d<-merge(bmi_exp_dat_clumped,c,by.x = "SNP",by.y = "MarkerName")
write.csv(d, file="outcome.csv") 

#将outcome数据读入TwoSampleMR
outcome_dat<-read_outcome_data(snps = bmi_exp_dat_clumped$SNP,
                               filename = "outcome.csv",
                               sep = ",",
                               snp_col = "SNP",
                               beta_col = "beta",
                               se_col = "se",
                               effect_allele_col = "effect_allele",
                               other_allele_col = "other_allele",
                               pval_col = "p")

#协同
dat<-harmonise_data(exposure_dat = bmi_exp_dat_clumped,
					outcome_dat = outcome_dat)
write.csv(dat, file="harmonized_data.csv") 

####3、进行MR#### 
mr(dat)
generate_odds_ratios(mr_res = mr(dat)) 


mr_method_list()#查看TwoSampleMR包里的mr方法
mr(dat,method_list = c("mr_ivw","mr_egger_regression","mr_weighted_median"))#选择所用的MR方法

#结果可视化
mr_scatter_plot(mr_results = mr(dat,method_list = c("mr_ivw","mr_egger_regression","mr_weighted_median")),dat)

#异质性检测
mr_heterogeneity(dat)

#run_mr_presso(dat,NbDistribution =3000)
#异质性可视化，对称，异质性越小
mr_funnel_plot(singlesnp_results = mr_singlesnp(dat))

#多效性检测，使用egger_intercept,评估截距是否大于0.05
mr_pleiotropy_test(dat)

mr_leaveoneout_plot(leaveoneout_results = mr_leaveoneout(dat))

Harmonization的作用：

1、将暴露和结局的SNP等位基因方向协同
2、根据EAF大小，剔除不能判断方向的palindromic SNP
3、剔除incompatible SNP（A/G vs. A/C）
协同之后需要看一下SNP与结局的P值，剔除p<5e-08的SNP。	剔除F值<10的SNP。

MR结果：

若结局是连续型变量（如：身高、体重），则可直接用beta表示；
若为二分类变量，则用OR表示。

异质性检测：mr_heterogeneity()
1、何为异质性

不同SNP之间的wald ratio 之间的差异，p<0.05，存在异质性

2、结果解读？
3、异质性有何影响
4、如果解决异质性？

run_mr_presso(dat,NbDistribution =3000)
结果p<0.05可能为离群值

统计效能：
结果好可放，不好不放

问题：
1、暴露数据和结局数据的SNO所对应的效应等位基因不一致（A/G vs. G/A）?

协同

2、一部分暴露-SNP在结局中找不到，怎么办？

答：看丢失数量；
丢失少可以不管，也可以找代理SNP（如100丢10几可接受，丢40几不可），
丢失多结果不可靠。

3、代理SNP（proxies）如何确定

答：同2，丢失多找代理SNP结果依然不可靠，
	丢失少是否找SNP可看个人习惯。
    找代理SNP网站：http://snipa.org/snipa3/
    找代理时，LD R²应尽可能大，一般设置为0.8。

结局数据选择注意要点：
1、是否有样本重叠(sample overlapping)？

样本最大重叠率在10%以下可接受？（不确定）
样本重叠率计算：
例：
		暴露和结局的GWAS_meta数据里面都包含了A队列（1千人），
		暴露的GWAS_meta数据有100万人，
		结局的GWAS_meta数据有50万人，
		则样本最大重叠率为： MAX(1千/100万,1千/50万 )。

2、结局数据中SNP量是否足够多？

SNP量最好在400万，500万以上