学习笔记-孟德尔随机化代码
1、暴露数据和结局数据的SNO所对应的效应等位基因不一致(A/G vs. G/A)?1、是否有样本重叠(sample overlapping)?2、一部分暴露-SNP在结局中找不到,怎么办?3、代理SNP(proxies)如何确定。2、结局数据中SNP量是否足够多?孟德尔随机化-----本地数据。4、如果解决异质性?
文章共1,666字 · 阅读需要大约6分钟
一键AI生成摘要,助你高效阅读
问答
·
孟德尔随机化-----本地数据
#打开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万以上
更多推荐
已为社区贡献1条内容
所有评论(0)