各位小伙伴们，什么样的MR分析才能入了审稿人的法眼呢？怎样才能从内卷中脱颖而出呢？想要自己的MR分析脱颖而出，当然更是卷工作量、分析方法和新颖度。Radial MR从可视化图入手，改变了原来散点图的横纵坐标，构建新的统计量，计算离群值的思路值得小伙伴们跟着果叔来学习一下。
双样本孟德尔随机化 （MR） 研究的汇总数据通常借助散点图进行可视化，其中单核苷酸多态性 （SNP） 结果关联与 SNP 暴露关联作图，以提供每个个体变体的因果效应估计的结果。将因果效应的标准逆方差加权 （IVW） 估计值叠加为拟合斜率也很方便，以查看单个 SNP 是否提供了支持或与总体共识相冲突的证据。不幸的是，传统的散点图并不是实现这一目标的最合适方法，因为SNP-结果关联以不同程度的精度进行估计。Radial-MR消除了对遗传数据进行重新编码的需要，并能够更直接地检测异常值和有影响力的数据点。Radial-MR提出了在进行 MR 研究时可以操作的更通用的建模框架，包括一种新形式的 MR-Egger 回归。
Radial-MR的优势：
1.提供双样本孟德尔随机化 （MR） 研究的汇总数据通常借助散点图进行可视化。散点图还用于解释标准逆方差加权 （IVW） 估计值和多效性鲁棒方法（如 MR-Egger 回归）的有效性。
2.径向图消除了对汇总数据进行预处理的需要（MR-Egger的先决条件），改进了IVW或MR-Egger分析中异常值和影响数据点的检测，并且可以合并用户所需的任何权重集。
3.提出了一种新形式的MR-Egger 回归形式
跑代码时卡顿、电脑不给力让人抓狂！找果叔试用稳定高速的服务器，让分析顺畅无比！
代码学不会？bug 频繁出现，束手无策？实操生信分析课程赶快学起来！滴滴果叔领取体验课程哦~
线上课程教学
课题设计、定制生信分析
云服务器租赁
加微信备注99领取使用
代码教程
加载R包
library(tidyverse)library(TwoSampleMR)library(gt)library(RadialMR)library(phenoscanner)library(reactable)
执行Radial MR
RadialMR 软件包要求的输入数据格式与 TwoSampleMR 不同，因此我们首先要将统一数据转换为TwoSampleMR::dat_to_RadialMR()TwoSampleMRRadialMR 所要求的格式。
# 格式化数据radial_dat <- mr_dat %>% filter(mr_keep == T) %>% dat_to_RadialMR()rmr_dat <- radial_dat$TC.AD# 运行 Radial MRbonff = 0.05/nrow(rmr_dat) # Bonferonni Correctionradial_ivw_res <- ivw_radial(rmr_dat, alpha = bonff)radial_egger_res <- egger_radial(rmr_dat, alpha = bonff)
与我们的核心MR方法一样，我们观察到径向MR因果估算表明，总胆固醇与 AD 风险增加有因果关系,存在显著的异质性。
通过检查径向图，我们还进一步观察到一些变异体被归类为异常值，这些异常值是根据其对全局异质性的贡献来确定的，用 Cochran 的 Q 统计量和 bonferonni 校正的显著性阈值进行量化。
可视化
ivw_radial_p <- plot_radial(radial_ivw_res, radial_scale = F, show_outliers = F)egger_radial_p <- plot_radial(radial_egger_res, radial_scale = F, show_outliers = F)cowplot::plot_grid(ivw_radial_p + coord_fixed(ratio=0.25) + theme(legend.position = 'bottom'),egger_radial_p + theme(legend.position = 'bottom'),align = 'h')
IVW Radial MR 方法确定了六个变异为潜在异常值。我们可以查询 PhenoScanner，以获得有关这些 SNP（以及与它们有 LD 关系的变异）的更多信息，包括它们与哪些其他性状相关。
去除异常值
phewas_dat <- phenoscanner(snpquery=radial_ivw_res$outliers$SNP,proxies='EUR',r2 = 0.8)radial_ivw_res$outliers %>%left_join(filter(phewas_dat$snps, snp == rsid), by = c("SNP" = "snp")) %>%select(SNP, Q_statistic, p.value, hg19_coordinates, a1, a2, hgnc) %>%reactable()
通过研究 PheWAS 的结果，我们可以发现三个异常值（rs7412、rs7568769 和rs8103315）都位于 APOE 区域内，并与阿尔茨海默病有关。这些变异与之前观察到的驱动 " Leave-one-out "分析中 IVW 因果估计的变异相同。
既然已经确定了异常值，我们就可以在去除异常值后重新运行核心 MR 方法。我们将Radial-MR结果加入统一数据集，并修改变量 mr_keep，以标记异常值进行清除。
## 修改 mrkeep 变量，将变体标记为径向 MR 的异常值，以便移除mr_dat_outlier <- mr_dat %>%left_join(radial_ivw_res$dat) %>%mutate(mr_keep = case_when(Outliers == "Outlier" ~ FALSE,TRUE ~ mr_keep))
剔除异常值后，IVW 估计值不再显著。
无异常值的 MR 分析
#标准化mr_res <- mr(mr_dat, method_list = c(c("mr_ivw_fe", "mr_egger_regression", "mr_weighted_median", "mr_weighted_mode")))## MR 分析mr_res_outlier <- mr(mr_dat_outlier, method_list = c("mr_ivw_fe", "mr_egger_regression", "mr_weighted_median", "mr_weighted_mode"))# 异质性统计mr_het_outlier <- mr_heterogeneity(mr_dat_outlier, method_list = c("mr_egger_regression", "mr_ivw"))mr_plt_outlier <- mr_pleiotropy_test(mr_dat_outlier)# Radial MRradial_dat_outlier <- mr_dat_outlier %>% filter(mr_keep == T) %>% dat_to_RadialMR()radial_res_outlier <- ivw_radial(radial_dat_outlier$TC.AD, alpha = 0.05/nrow(radial_dat_outlier$TC.AD))bind_rows(mr_res %>% mutate(model = "w/ Outliers"),mr_res_outlier %>% mutate(model = "w/o Outliers")) %>%select(model, method, nsnp, b, se, pval) %>%gt(groupname_col = "model") %>%fmt_number(columns = c("b", "se")) %>%fmt_number(columns = pval,rows = pval > 0.001,decimals = 3) %>%fmt_scientific(columns = pval,rows = pval <= 0.001,decimals = 1)# 可视化scatter_p <- mr_scatter_plot(mr_res, mr_dat)scatter_out_p <- scatter_p[[1]] + theme_bw() +guides(color=guide_legend(nrow =1)) +labs(title = "w/ outliers") +theme(text = element_text(size = 8),legend.position = 'bottom')scatter_outlier_p <- mr_scatter_plot(mr_res_outlier, mr_dat_outlier)scatter_outlier_out_p <- scatter_outlier_p[[1]] + theme_bw() +guides(color=guide_legend(nrow =1)) +labs(title = "w/o outliers") +theme(text = element_text(size = 8),legend.position = 'bottom')radial_outlier_p <- plot_radial(radial_res_outlier, radial_scale = F, show_outliers = F)cowplot::plot_grid(scatter_out_p + theme(legend.position = 'none'),scatter_outlier_out_p + theme(legend.position = 'none'),ivw_radial_p + coord_fixed(ratio=0.25) + theme(legend.position = 'none') + labs(title = "w/ outliers"),radial_outlier_p + coord_fixed(ratio=0.2) + theme(legend.position = 'none') + labs(title = "w/o outliers"))
APOE 区域
在初步分析中，Radial-MR 确定了六个 SNP 为潜在异常值。其中三个 SNP 位于 APOE 区域，该区域以其效应量大和多效应性而著称，根据其 Q 值统计，这些 SNP 显然是异常值。相比之下，其他三个变异位于基因组的其他地方，其 Q 值明显较小。个别 APOE 变体对全局异质性统计量的巨大贡献可能会使其他变体的 Q 统计量产生偏差。
无 APOE 的Radial MR
mr_dat_wo_apoe <- mr_dat %>%mutate(apoe_region = case_when(chr.outcome == 19 & between(pos.outcome, 44912079, 45912079) ~ TRUE,TRUE ~ FALSE),gws.outcome = ifelse(pval.outcome < 5e-8, TRUE, FALSE),mr_keep = case_when(apoe_region == TRUE ~ FALSE,gws.outcome == TRUE ~ FALSE,TRUE ~ mr_keep))# Radial MRradial_dat_wo_apoe <- mr_dat_wo_apoe %>% filter(mr_keep == T) %>% dat_to_RadialMR()radial_res_wo_apoe <- ivw_radial(radial_dat_wo_apoe$TC.AD, alpha = 0.05/nrow(radial_dat_wo_apoe$TC.AD))radial_wo_apoe_p <- plot_radial(radial_res_wo_apoe, radial_scale = F, show_outliers = F)
在进行Radial MR分析前移除 APOE 区域后，只有两个 SNP 被标记为异常值。这两个SNP分别是rs6504872和rs9391858，rs6504872在移除APOE区域后仍被认为是异常值，而rs9391858则在这一步骤后成为异常值。这表明，包含 APOE 区域会使其他变异体的 Q 统计量出现偏差。
可视化
cowplot::plot_grid(ivw_radial_p + coord_fixed(ratio=0.25) + theme(legend.position = 'none') + labs(title = "w/ APOE"),radial_wo_apoe_p + coord_fixed(ratio=0.2) + theme(legend.position = 'none') + labs(title = "w/o APOE"),align = 'h')
小结
Radial-MR 改进了MR分析的可视化解释，并显示对 Cochran's Q 统计量贡献较大的数据点，这些数据点很可能是异常值。根据这一统计标准，可以从下游MR分析中剔除个别变异体，以获得更可靠的因果关系。不过，在采取剔除所有异常值的策略时应小心谨慎，直到几乎不存在异质性。这一点在考虑 APOE 区域时尤为突出，尽管该区域的 Q 统计量适中，但仍有其他变异体被归类为异常值。排除特定 SNP 时，还应考虑该 SNP 是否与代表结果多效应途径的单独表型相关。最后，大家如果对生信分析感兴趣但还不熟悉，又想尝试一下处理自己的数据，不妨试一下果叔开发的生信云平台哦，一键出图，一键导出CNS级别的Figture！！赶快去试试吧！！点击 http://www.biocloudservice.com/home.html。
参考文献
1.Bowden, J., et al., Improving the visualization, interpretation and analysis of two-sample summary data Mendelian randomization via the Radial plot and Radial regression. International Journal of Epidemiology, 2018. 47(4): p. 1264-1278.
2.https://andrewslabucsf.github.io/MR-tutorial/scripts/radial_mr.html
不会分析还想用生信工具助力发文咋办？有这顾虑的朋友，想一步到位就带着想法来，不论是代码实操还是在线文章结果复现，果叔照样能提供，还有大家都想要的服务器，找果叔获取就对了！
往期回顾
01 |
02 刚逃过“三花淡奶”又陷入“预制菜”！“NHANES最新数据+多变量回归”实锤：超加工食品会加重骨质疏松！不爱做饭的亲，注意咯！ |
03 |
04 思路简直开挂，已经复现坐等毕业了！青岛大学的网络药理学研究咋就这么牛？看来“药” 做就做高级范，毕业才能so easy！ |
