孟德尔随机化：代码分享（二）

生信菜鸟团

发布于 2023-09-09 17:03:01

3.2K00

代码可运行

文章被收录于专栏：生信菜鸟团生信菜鸟团

运行总次数：0

代码可运行

第一部分在这里：孟德尔中介分析全流程代码（一）

话不多说，上代码——

1加载包

library(tidyverse)
library(data.table)
library(ivpack)
library(meta)
library(devtools)
library(pacman)
library(MendelianRandomization)
library(MRInstruments)
library(ieugwasr)
library(phenoscanner)
library(LDlinkR)
library(mr.raps)
library(MRPRESSO)
library(extrafont)
library(anchors)

2自定义一个函数

留心看——load(paste0(deparse(substitute(exposure)), "_", deparse(substitute(outcome)), ".RData")) 这部分的数据就是上一次运行的结果，记得命名要规范，是以exposure_outcome.RData格式命名的，根据自己的暴露和结局修改一下即可。

#'Creates table of IVW results for the exposure-outcome combination
ivw_table = function(exposure, outcome) {
  load(paste0(deparse(substitute(exposure)), "_", deparse(substitute(outcome)), ".RData"))
  dat <- dat[dat$mr_keep == TRUE, ]##只有mr_keep是TRUE的SNP才是真正用于计算MR结果的IV
  mr_object = mr_input(bx = dat$beta.exposure, 
                       bxse = dat$se.exposure, 
                       by = dat$beta.outcome, 
                       byse = dat$se.outcome, 
                       exposure = deparse(substitute(exp)), 
                       outcome = deparse(substitute(out)), 
                       snps = dat$SNP)
  ivw_res = mr_ivw(mr_object) ##熟悉的函数对不对

mr_ivw函数

mr_ivw 函数实现了反方差法，非正式地称为 "托比-约翰逊 "法。对于单一遗传变异，这就是简单的比率法。ivw法也是MR分析的核心哦~

mr_ivw(
  object,
  model = "default",
  robust = FALSE,
  penalized = FALSE,
  weights = "simple",
  psi = 0,
  correl = FALSE,
  distribution = "normal",
  alpha = 0.05,
  ...
)

object MRInput 对象。 model 应使用的模型类型："默认"、"随机 "或 "固定"。随机效应模型（"random"）是一种乘法随机效应模型，允许加权线性回归中的过度分散（残差标准误差不固定为 1，但不允许取低于 1 的值）。固定效应模型（"fixed"）将残差标准误差设为 1。默认 "设置是在有 3 个或更少的遗传变异时使用固定效应模型，否则使用随机效应模型。 robust 表示该方法是否应使用 robustbase 软件包中的 lmrob() 函数进行稳健回归，而不是使用标准线性回归（lm）。 penalized 表示是否应对权重进行惩罚，以降低具有离散比率估计值的遗传变异对分析的贡献。 weights 在加权回归中使用的权重。如果使用 "simple"（简单）（默认选项），则 IVW 估计值等同于根据比率估计值方差的最简单表达式（delta 扩展的一阶项--与结果相关性的标准误差除以与暴露相关性），使用逆方差权重对每个变异的比率估计值进行元分析。如果是 "delta"，则方差表达式为 delta 扩展的二阶项。二阶项包含了遗传与暴露相关性的不确定性--使用简单加权法可以忽略这种不确定性。 psi 基因与暴露的相关性与样本重叠产生的每个变异体与结果的相关性之间的相关性。默认值为 0，相当于严格的双样本孟德尔随机分析（无重叠）。如果样本之间完全重叠，则相关性应设置为暴露与结果之间的观察相关性。只有在权重选项设置为 "delta "时，该相关性才会用于计算标准误差。 correl 如果基因变异是相关的，则可以考虑这种相关性。必须在 MRInput 对象中提供相关性矩阵：该矩阵的元素是各个变异体之间的相关性（对角元素为 1）。如果 MRInput 对象中指定了相关性矩阵，则 correl 设置为 "true"。如果 correl 设置为 "true"，则稳健值和惩罚值均为 "false"，权重设置为 "simple"。 distribution 用于计算置信区间的分布类型。可选项有 "正态分布"（默认）或 "t-分布"。 alpha 用于计算置信区间的显著性水平。默认值为 0.05。

接下来就是置信区间和p值的计算啦~

 if (deparse(substitute(outcome) == "LOAD")) {##这里也要换成自己的结局！
    ivw_res@Estimate = exp(ivw_res@Estimate) ##exp()函数用于计算e的幂，即e^y或者我们可以说y的 index
    ivw_res@CILower = exp(ivw_res@CILower)
    ivw_res@CIUpper = exp(ivw_res@CIUpper)
  }
  
  if (deparse(substitute(exposure) == "LOAD")) {##这里也要换成自己的结局！
    ivw_res@Estimate = log(2)*(ivw_res@Estimate)
    ivw_res@CILower = log(2)*(ivw_res@CILower)
    ivw_res@CIUpper = log(2)*(ivw_res@CIUpper)
  }
  
  if (ivw_res@Pvalue < 0.001) {
    df = data.frame(exposure = deparse(substitute(exposure)), 
                    outcome = deparse(substitute(outcome)), 
                    N_SNPs = ivw_res@SNPs,
                    Estimate = paste0(sprintf("%.2f", round(ivw_res@Estimate, 2)), " (", sprintf("%.2f", round(ivw_res@CILower, 2)), ", ", sprintf("%.2f", round(ivw_res@CIUpper, 2)), ")"),
                    P = format(ivw_res@Pvalue, digits = 3, scientific = TRUE)
    )
  } else {
    df = data.frame(exposure = deparse(substitute(exposure)), 
                    outcome = deparse(substitute(outcome)), 
                    N_SNPs = ivw_res@SNPs,
                    Estimate = paste0(sprintf("%.2f", round(ivw_res@Estimate, 2)), " (", sprintf("%.2f", round(ivw_res@CILower, 2)), ", ", sprintf("%.2f", round(ivw_res@CIUpper, 2)), ")"),
                    P = sprintf("%.3f", round(ivw_res@Pvalue, 3))
    )
  }
  return(df)
}

3运行这个函数

yourdat= ivw_table(yourexposure, youroutcome)

如果你的结局是由多个表型或性状构成的，那就这样：

EA_IDPs_ivw = rbind(
  ivw_table(EA, SA),
  ivw_table(EA, vol),
  ivw_table(EA, CT),
  ivw_table(EA, LGI),
  ivw_table(EA, MC),
  ivw_table(EA, IC),
  ivw_table(EA, FA_cortical),
  ivw_table(EA, MD_cortical),
  ivw_table(EA, ICVF),
  ivw_table(EA, ODI),
  ivw_table(EA, FA_WM),
  ivw_table(EA, MD_WM),
  ivw_table(EA, hippocampus_vol_left),
  ivw_table(EA, hippocampus_vol_right),
  ivw_table(EA, WMH_vol)
)

加一下标签：

EA_IDPs_ivw$outcome = factor(
  x = EA_IDPs_ivw$outcome,
  levels = c(
    "SA", 
    "vol", 
    "CT", 
    "LGI", 
    "MC", 
    "IC", 
    "FA_cortical", 
    "MD_cortical", 
    "ICVF", 
    "ODI", 
    "FA_WM", 
    "MD_WM", 
    "hippocampus_vol_left",
    "hippocampus_vol_right", 
    "WMH_vol"
    ), 
  labels = c(
    "Surface area", 
    "Volume",
    "Cortical thickness", 
    "Local gyrification index", 
    "Mean curvature", 
    "Intrinsic curvature", 
    "Fractional anisotropy (cortical)", 
    "Mean diffusivity (cortical)", 
    "Intracellular volume fraction", 
    "Orientation dispersion index", 
    "Fractional anisotropy (white matter tracts)", 
    "Mean diffusivity (white matter tracts)", 
    "Volume of left hippocampus", 
    "Volume of right hippocampus", 
    "White matter hyperintensities volume"
  )
)

后期加上作图函数以后，就会得到下面这张森林图：