TwoSampleMR-R教程 两样本孟德尔随机化(原来真的就是这么简单……)
1. 工具变量
install.packages("remotes")
remotes::install_github("MRCIEU/TwoSampleMR")##To update the package just run the command again
install.packages("readxl")
读取本地文件
- 方法1:先改变量名,然后
format_data()
读取自己整理好的工具变量文件,命名如下图。注意,命名区分大小写,不正确则无法读取
##name:SNP beta se effect_allele eaf exposure chr position gene pval samplesize mr_keep
HcyIV <- as.data.frame(read_excel("IV_Hcy.xlsx"))
library(TwoSampleMR)
Hcy_exp_dat <- format_data(HcyIV, type="exposure")
- 方法2:直接用
read_exposure_data()
读,然后code改。尽可能提供多的信息,方便后期画图
HcyIV.2 <- read_exposure_data(filename = "IV_Hcy2.csv",sep = ",",snp_col = "SNP",beta_col = "beta",se_col = "se",effect_allele_col = "effect_allele",#other_allele_col = "a2",eaf_col = "eaf",pval_col = "pval",#units_col = "Units",#gene_col = "Gene",samplesize_col = "samplesize"
)
- 最后改好自己的暴露名
Hcy_exp_dat$exposure <- "Homocysteine"
工具变量必须为独立变量。根据参考文献选择ref人群和r2
## Clumping via IEU GWAS database
Hcy_exp_dat_clumped <- clump_data(Hcy_exp_dat, clump_r2 = 0.2, pop="EUR")
# default using EUR population from 1000 Genomes project
## Read reference for suitable r2
2. 结局变量
选好工具变量之后,就要将工具变量从结局的GWAS中提取出来。
可以像上面一样从本地读取,也可以从R包里提。下面展示从R包提取的方法。
从数据库获取结局相关的GWAS
## Get GWAS data from the IEU GWAS database that rooted in the package
ao <- available_outcomes()
#> API: public: http://gwas-api.mrcieu.ac.uk/
ao_outcome <- ao[grepl("stroke", ao$trait), ]
View(ao_outcome)
一个trait很多实验,一个实验一个id;实验来自不同人群和biobank,选择最新的样本量最大的高质量的研究(看year and author)、选择相应研究的id。
选择Ischemic stroke的id是
ieu-a-1108
也可以直接到网站查找,可以顺便查看id Batch的细节。
https://gwas.mrcieu.ac.uk/
但更建议手动查找,因为网站不全!
从结局GWAS提取工具变量
IS_out_dat <- extract_outcome_data(snps = Hcy_exp_dat_clumped$SNP, # use data after clumpingoutcomes = 'ieu-a-1108',proxies = TRUE, #defaultrsq=0.8, #default; min is 0.6palindromes = 1, #Allow palindromic SNPs? Default is 1 (yes)maf_threshold = 0.3 #If palindromes allowed then what is the maximum minor allele frequency of palindromes allowed? Default is 0.3.
)
软件会自动根据LD选取proxy,default R2=0.8
Extracting data for 18 SNP(s) from 1 GWAS(s)
Finding proxies for 1 SNPs in outcome ieu-a-1108
Extracting data for 1 SNP(s) from 1 GWAS(s)
注意网上extract默认允许palidromes。Palidrome是啥?
A palindromic SNP (also known as an ‘ambiguous’ SNP) is a SNP in which the possible alleles for that SNP pair with each other in the double helix structure
Palidrome就是“回文序列”,简单来说就是,A/T、C/G这样的结构。如图所示,A1=A,A2=T,无法根据A1和A2来判断暴露和结局的两个SNP是否方向一致(因为A对应就是T,T对应就是A)。此时就需要根据allele frequency来判断。如下图,暴露的eaf=0.11,如果结局也是和暴露一个方向,其eaf也应该在0.11左右(默认可接受差距=0.3);但是其eaf=0.91,太大了。所以我们推断,结局allele的方向和暴露的相反,所以其效应effect应该反过来= - (-0.05)=0.05
注意:如结局summary是本地文件,建议根据IV,merge
得到结局summary,然后手动检查和判断回文结构。因为下一步harmonize是没有再allow palindromic的选项了。
harmonising的时候默认所有序列都是顺读的(presented on the forward strand)。
3. 数据预处理
对数据进行预处理,使其效应等位与效应量保持统一,就是调整exp和out上位点的方向和beta,使其统一。
## To harmonise the
dat <- harmonise_data(exposure_dat = Hcy_exp_dat, outcome_dat = IS_out_dat)
4. MR和敏感性分析等
## MR
res <- mr(dat)## Heterogeneity statistics
mr_heterogeneity(dat)# pleiotropy test
mr_pleiotropy_test(dat)# leave-one-out analysis
res_loo <- mr_leaveoneout(dat)
默认五种方法:MR Egger、Weighted median、Inverse variance weighted、Simple mode和Weighted mode
用mr_method_list()
查看当前支持的其他统计方法。
在mr函数中添加想要使用的方法
mr(dat, method_list=c("mr_egger_regression", "mr_ivw"))
5. 可视化
# scatter plot
p1 <-mr_scatter_plot(res, dat)
length(p1) # to see how many plots are there
p1[[1]]# forest plot
res_single <- mr_singlesnp(dat)#,all_method=c("mr_ivw", "mr_two_sample_ml")) to specify method used
p2 <- mr_forest_plot(res_single)
p2[[1]]## funnel plot for all
p3 <- mr_funnel_plot(res_single)
p3[[1]]
以上这些函数都可以同通过method=
增加特定的方法。
导出图片
### save images ##########
library(ggplot2)
ggsave(p1[[1]], file="mr_scatter_plot.png", width = 7, height=7, dpi=900)
6. 结果整合
一键生成所有MR分析、敏感性分析和画图,并保存在html、word或pdf。默认在工作目录生成一个html。
mr_report(dat)
但是我没成功,不推荐,还是每个分析和图都手动做吧。一键生成不利于发现问题。
还可以用这个方法:
## To combine all resultsres<-mr(dat)
het<-mr_heterogeneity(dat)
plt<-mr_pleiotropy_test(dat)
sin<-mr_singlesnp(dat)
all_res<-combine_all_mrresults(res,het,plt,sin,ao_slc=T,Exp=T,split.exposure=F,split.outcome=T)
head(all_res[,c("Method","outcome","exposure","nsnp","b","se","pval","intercept","intercept_se","intercept_pval","Q","Q_df","Q_pval","consortium","ncase","ncontrol","pmid","population")])
7. Variance explained
out <- directionality_test(dat)
kable(out)
MR假定工具变量先影响暴露,然后通过暴露影响结局,但也有可能是一个cis-acting SNP先影响了基因表达或DNA甲基化。所以,为了确定这个假定的方向性,我们使用 the Steiger test
分别计算工具变量对暴露和结局的variance explained,并判断结局的variance是否小于暴露;如果是,则方向正确。
这个测试有可能因为研究GWAS本身的质量不够好而导致结果不可靠,可以再做sensitivity analyses。方法如下:
- Provide estimates of measurement error for the exposure and the outcome, and obtain an adjusted estimate of the causal direction.
- For all possible values of measurement error, identify the proportion of the parameter space which supports the inferred causal direction
mr_steiger(p_exp = dat$pval.exposure, p_out = dat$pval.outcome, n_exp = dat$samplesize.exposure, n_out = dat$samplesize.outcome, r_xxo = 1, r_yyo = 1,r_exp=0)
Citation
TwoSampleMR的R包: Hemani G, Zheng J, Elsworth B, et al. The MR-Base platform supports systematic causal inference across the human phenome. Elife 2018;7
directionality test: Hemani G, Tilling K, Davey Smith G. Orienting the causal relationship between imprecisely measured traits using GWAS summary data. PLoS Genet. 2017 Nov 17;13(11):e1007081. doi: 10.1371/journal.pgen.1007081. Erratum in: PLoS Genet. 2017 Dec 29;13(12 ):e1007149. PMID: 29149188; PMCID: PMC5711033.
有问题请直接评论,不回复私信。
TwoSampleMR-R教程 两样本孟德尔随机化(原来真的就是这么简单……)相关推荐
- R数据分析:孟德尔随机化实操
好多同学询问孟德尔随机化的问题,我再来尝试着梳理一遍,希望对大家有所帮助,首先看下图1分钟,盯着看将下图印在脑海中: 上图是工具变量(不知道工具变量请翻之前的文章)的模式图,明确一个点:我们做孟德尔的 ...
- R数据分析:孟德尔随机化分析文献解析和实例操练
最近抽空研读了一篇探讨高血压和肾功能关系的文献,记录下来分享给大家,主要也是想看看孟德尔随机化的统计分析结果在论文中是如何呈现的,之后我会给大家写写孟德尔随机化的统计分析在R语言中的做法,希望可以帮助 ...
- spss和sas和python_T检验第三篇(SPSS,SAS,R,Python) 两样本T检验
两样本T检验,和终于来到T检验的最后一个章节,两样本T检验. 两样本T检验的应用条件为:1.独立的随机样本 2.资料应当服从正态分布 3.方差齐性 即我们要在前面两种T检验的前提下,做多一个方差齐性检 ...
- R语言丨根据VCF文件设计引物,自动识别两样本差异SNP位点,调用samtools获取上下游参考序列
根据变异位点设计引物序列 今天碰到一个新问题:假如有一个vcf文件储存了两个样品的变异位点基因型数据,每行代表一个位点,我现在想找出两样本差异的SNP位点,再把差异位点用[REF/ALT]的形式表示, ...
- [文献分享] 父母炎症性肠病与儿童自闭症(国家登记数据队列研究、连锁不平衡分数回归、多基因风险评分、孟德尔随机化)
文献来源:Sadik A, Dardani C, Pagoni P, et al. Parental inflammatory bowel disease and autism in children ...
- 一起来学孟德尔随机化(Mendelian Randomization)
孟德尔随机化最近实在是太火了,想不关注都不行,最近也花了点时间研究了一下,和大家分享一下,共同学习. 什么是孟德尔随机化? 在19世纪,孟德尔用豌豆花作为实验材料,通过对豌豆花颜色.形状等特征的观察和 ...
- 解析基因影响:孟德尔随机化的创新思维
一.引言 在当今的遗传学和生物学研究中,我们对基因对个体特征和性状的影响的理解变得更加深入.然而,基因影响的复杂性和多样性给我们带来了巨大的挑战.为了更好地揭示基因影响的本质和机制,我们需要采用创新的 ...
- 非参数统计:两样本和多样本的Brown-Mood中位数检验;Wilcoxon(Mann-Whitney)秩和检验及有关置信区间;Kruskal-Wallis秩和检验
目录 两样本和多样本的Brown-Mood中位数检验 例3.1我国两个地区一些(分别为17个和15个)城镇职工的工资(元): Wilcoxon(Mann-Whitney)秩和检验及有关置信区间 例3. ...
- 孟德尔随机化分析时,异质性太强怎么办
孟德尔随机化分析是用来研究不同个体之间的差异的统计方法,如果样本中的个体太异质,可能会导致结果的不确定性增加.在这种情况下,可以考虑以下几种方法来减少异质性的影响: 对样本进行分层,即按照相关特征将样 ...
最新文章
- NanoPi NEO Air使用二:固件烧录
- mysql的中文乱码url,MySQL 中文显示乱码
- CLion 输出遇到乱码解决办法,GBK和utf-8的转换
- 浙大 PAT b1029
- 浙大通讯与计算机网络离线作业,2015浙大通信与计算机网络离线作业.doc
- python 基础及资料汇总
- 一种多功能语音识别技术和音乐播放器相结合的方法
- zabbix 服务器监控之数据库操作
- 计算机创新设计2大赛获奖作品3Done,走向3D创意世界——3Done创客设计比赛
- 用matlab做仿真实验难不难,SIMULINK仿真实验心得体会
- java根据word模板生成word文档_根据Word模板生成Word文件 (JAVA POI)
- 自己拥有一台服务器可以做哪些很酷的事情
- Apple中文社区平台的 Mac 用户群体
- shell输出标准时间格式
- 【原创】VBA学习笔记(15)VBA的参数传递:ByVal 和 ByRef 的区别
- 美味连连-QQ游戏辅助-简单实用的QQ游戏美味连连辅助(非外挂)
- iOS 分析一次有意思的需求——HTML代码注入
- RTA PAVIA CSD ET04控制器
- 超级牛批的IP地址查询工具
- MySQL 和 Oracle 大数据量分页查询方法及其优化
热门文章
- 2020 China Collegiate Programming Contest, Weihai B Labyrinth
- springboot模板整合(四)邮箱验证
- 程序员真是太太太太太有趣了!!!
- html5 video 停止播放视频,html5 video怎么停止播放视频
- Android 开发 之 8.0应用安装权限(未知应用权限安装)
- Yield Guild Games: RON 质押来啦!
- 小程序再升级:支持创建微信小店小程序
- 20+应急响应工具包合集(附下载地址)
- 计算机图形学试题题库,计算机图形学题库及答案
- Nacos的动态配置源码解析