使用TwoSampleMR进行孟德尔随机化分析

原创

凑齐六个字吧

发布于 2026-02-24 09:44:38

1520

TwoSampleMR 是一个用于进行双样本孟德尔随机化分析的R包，由 MRC Integrative Epidemiology Unit开发，并与OpenGWAS数据平台深度整合。主要用于因果推断分析，来评估某个“暴露因素”是否影响某个“结局”。

分析流程

本研究使用的暴露与结局 GWAS 汇总数据均来自OpenGWAS 数据库，其中暴露因素为甘油三酯（数据库 ID：ieu‑b‑111），结局为结直肠癌（数据库 ID：ebi‑a‑GCST90018808）。

1.导入

rm(list=ls())
library(TwoSampleMR)
library(VariantAnnotation) 
library(gwasglue) 
library(dplyr)
library(tidyr)
library(CMplot)

2.暴露数据处理

进行孟德尔随机化分析时，需要提供一个包含暴露工具变量的数据框。数据框中的每一行代表某一个暴露的一个遗传变异的相关信息。

进行MR分析所需的最低限度信息包括以下内容：1.SNP：rs号（如 rs123456）；2. beta：效应值。若该性状为二分类变量（如患病/未患病），则应使用对数优势比（log(OR)）；3. se：效应值的标准误；4. effect_allele：与beta所示效应对应的等位基因。

以下信息虽非必需，但对MR分析非常有用，建议尽可能提供：1. other_allele：非效应等位基因；2. eaf：效应等位基因频率；3. Phenotype：该SNP所关联的表型（暴露）名称；

还可以提供以下额外信息：1. chr：变异所在的染色体编号；2. position：变异在染色体上的物理位置（碱基对位置）；3. samplesize：用于估计该效应值的总样本量；4. ncase：病例数（适用于二分类性状）；5. ncontrol：对照数（适用于二分类性状）；6. pval：该SNP与暴露关联的P值；7. units：效应值所使用的单位（如 mmol/L、kg 等）；8. gene：该SNP所在的基因或其他功能注释信息

1.把公共数据库中的GWAS数据转化成适用于TwoSampleMR分析的数据。

通常来说，从GWAS相关的数据库(比如openGWAS)中得到的数据可以通过函数转换转换成适合于TwoSampleMR分析的格式。

inputFile="ieu-b-111.vcf.gz"     
# 读取输入文件并
vcfRT <- readVcf(inputFile)
# 对输入文件进行格式转换
data=gwasvcf_to_TwoSampleMR(vcf=vcfRT, type="exposure")
head(data)

2.对数据按照pvalue进行过滤

# 根据pvalue<5e-08进行数据过滤
outTab <- subset(data, pval.exposure<5e-08)
write.csv(outTab, file="exposure.pvalure_filter.csv", row.names=F)

3.去除连锁不平衡的SNP

exposureFile="exposure.pvalure_filter.csv"     
# 读取数据
exposure_dat <- read_exposure_data(filename=exposureFile,
                       sep = ",",
                       snp_col = "SNP",
                       beta_col = "beta.exposure",
                       se_col = "se.exposure",
                       effect_allele_col = "effect_allele.exposure",
                       other_allele_col = "other_allele.exposure",
                       eaf_col = "eaf.exposure",
                       samplesize_col = "samplesize.exposure",
                       clump = F)

# 去除连锁不平衡的SNP
# 这里需要获取一下token
Sys.setenv(OPENGWAS_JWT="换成自己token")
## 连锁不平衡(LD)聚类
exposure_clumped <- clump_data(exposure_dat, 
                                   clump_kb=10000, # 默认10000kb范围内
                                   clump_r2=0.001, # r2阈值为0.001
                                   clump_p1 = 1, # 聚类显著性水平=1
                                   pop = "EUR"
                                   )
write.csv(exposure_clumped, file ="exposure_LD.csv", row.names=F)

4.去除弱工具变量

# 前面按照P值小于5*10e-8已经去除过一次
Ffilter=10        #F值过滤条件
    
#读取输入文件
inputFile="exposure_LD.csv"  
dat<-read.csv(inputFile, header=T, sep=",", check.names=F)

#计算F检验值
N=dat[1,"samplesize.exposure"]     #获取样品的数目
dat=transform(dat,R2=2*((beta.exposure)^2)*eaf.exposure*(1-eaf.exposure))     #计算R2
dat=transform(dat,F=(N-2)*R2/(1-R2))      #计算F检验值

#根据F值>10进行过滤, 删除弱工具变量
outTab=dat[dat$F>Ffilter,]
write.csv(dat, "exposure_F.csv", row.names=F)

5.去除混杂因素

LDlinkR方法经常失败

# 如果得到的SNP很多，就需要去除混杂因素相关的SNP；如果得到的SNP很少的话，就不需要去除混杂因素相关的SNP

## 1.LDlinkR方法## 这种方式经常失败
options(timeout = 6000000)
library(LDlinkR)     
inputFile="exposure_F.csv"          

# 1.读取输入文件
dat=read.csv(inputFile, header=T, sep=",", check.names=F)

# 2.SNP分组(LDlinkR要求每次分析小于50个SNP)
snpId=dat$SNP
y=seq_along(snpId)
chunks <- split(snpId, ceiling(y/50))

# 3.每次分析一个分组中的SNP
# 如果VPN影响了LDlinkR的链接，就需要去除
outTab=data.frame()
for(i in names(chunks)){
    # 混杂因素分析
    confounder=LDtrait(
      snps = chunks[[i]], # snp范围是1-50
        pop = "EUR",
        r2d = "r2",
        r2d_threshold = 0.1, # 一般设置默认数值即可
        win_size = 5e+05,
        token = "fd84165c4595",
        file = FALSE,
        genome_build = "grch37",
)
    outTab=rbind(outTab, confounder$results)
    Sys.sleep(1)  # 单位：秒
}
# 3.输出SNP相关形状的表格
write.csv(outTab, "confounder.result.csv", row.names=F)

## 2.FastTraitR方法###########
inputFile="exposure_F.csv"          
## 读取输入文件
dat_input=read.csv(inputFile, header=T, sep=",", check.names=F)
library(FastTraitR)
snp <- dat_input$SNP
write.csv(snp,file = "snp.txt")
res = look_trait(rsids=snp, pval=1e-5)

## 3.去除混杂因素
## 如果相应的snp出现在表型中那就需要去除
## 找到与脂质代谢相关的trait
keywords <- c(
  "adiponectin",
  "lipid",
  "triglyceride",
  "cholesterol",
  "LDL"
)
# 构建正则表达式（匹配任意一个关键词）
pattern <- paste(keywords, collapse = "|")
Snp_Del <- res$rsid[grepl(pattern, res$trait, ignore.case = TRUE)]
Snp_Del
#Snp_Del=c("rs13078960", "rs2030323")
dat=dat_input[!dat_input$SNP %in% Snp_Del,]
write.csv(dat, "exposure_confounder.csv", row.names=F)

6.结局数据处理

Name_exposure="Triglyceride"                  #图形中展示暴露数据的名称
Name_outcome="colorectal cancer"       #图形中展示结局数据的名称

exposureFile="exposure_confounder.csv"       #暴露数据输入文件
outcomeFile="ebi-a-GCST90018808.vcf.gz"              #结局数据输入文件

# 2.数据读取
exposure_dat<-read_exposure_data(filename=exposureFile,
                       sep = ",",
                       snp_col = "SNP",
                       beta_col = "beta.exposure",
                       se_col = "se.exposure",
                       effect_allele_col = "effect_allele.exposure",
                       other_allele_col = "other_allele.exposure",
                       eaf_col = "eaf.exposure",
                       clump = F)

# 3.读取结局数据的vcf文件，并对数据进行格式转换
vcfRT <- readVcf(outcomeFile)
outcomeData=gwasvcf_to_TwoSampleMR(vcf=vcfRT, type="outcome")

# 4.从结局数据中提取工具变量，只保留在两个数据框中都存在的SNP
outcomeTab<-merge(exposure_dat, outcomeData, 
                  by.x="SNP", by.y="SNP")
## 保留跟结局有关的列
write.csv(outcomeTab[,-(2:ncol(exposure_dat))], file="outcome.csv")

# 5.读取整理好的结局数据
outcome_dat<-read_outcome_data(snps=exposure_dat$SNP,
                 filename="outcome.csv", sep = ",",
                 snp_col = "SNP",
                 beta_col = "beta.outcome",
                 se_col = "se.outcome",
                 effect_allele_col = "effect_allele.outcome",
                 other_allele_col = "other_allele.outcome",
                 pval_col = "pval.outcome",
                 eaf_col = "eaf.outcome")

# 6.将暴露数据和结局数据合并
exposure_dat$exposure=Name_exposure
outcome_dat$outcome=Name_outcome
dat<-harmonise_data(exposure_dat=exposure_dat,
                    outcome_dat=outcome_dat)

# 7.输出用于孟德尔随机化的工具变量
dat[dat$mr_keep=="TRUE",] # 找到是否真正符合MR的数据
write.csv(outTab, file="table.SNP.csv", row.names=F)

# MR-PRESSO异常值检测
## 找到偏倚的SNP
## 如果Pvalue小于0.05就是偏倚的SNP
presso <- run_mr_presso(dat)
write.csv(presso[[1]]$`MR-PRESSO results`$`Outlier Test`, file="PRESSO.csv")

# 9.孟德尔随机化t=mr(dat)
mrResult=mr(dat)

#  mr_method_list()$obj # 选择孟德尔随机化的方法
# mrResult=mr(dat, method_list=c("mr_ivw", "mr_egger_regression", "mr_weighted_median", "mr_simple_mode", "mr_weighted_mode"))

# 对结果进行OR值的计算
mrTab <- generate_odds_ratios(mrResult)
write.csv(mrTab, file="MRresult.csv", row.names=F)

# 10.异质性分析
heterTab=mr_heterogeneity(dat)
write.csv(heterTab, file="heterogeneity.csv", row.names=F)

# 11.多效性检验
pleioTab=mr_pleiotropy_test(dat)
write.csv(pleioTab, file="pleiotropy.csv", row.names=F)

7.可视化

# 12.绘制散点图
pdf(file="scatter_plot.pdf", width=7.5, height=7)
mr_scatter_plot(mrResult, dat)
dev.off()

横轴：每个SNP对暴露因素的效应大小。纵轴：同一SNP对结局的效应大小。里面的每一个点代表了snp，十字图形代表了效应值的范围，横线代表了snp对暴露因素的影响范围，竖线代表了snp对结局的影响范围。然后其中五种分析方法得到了5条线段，如果线段呈左下到右上的趋势则说明甘油三酯能够增加结直肠癌的发病风险，如果线段如图所示则说明两者并没有明显的影响关系。

# 13.森林图
res_single=mr_singlesnp(dat)      
pdf(file="forest.pdf", width=7, height=6.5)
mr_forest_plot(res_single)
dev.off()

纵坐标代表工具变量snp，横坐标代表着效应值。每个snp右边对应的横线和点代表了效应值和效应值的波动范围。在0值左边的snp代表是保护因素，在0值右边的snp代表是危险因素。其中红色的是展示了孟德尔随机化下分析的综合效应，在这里两种分析方法的总和效应都横跨过了0，说明暴露因素与结局并不相关。

# 14.漏斗图
pdf(file="funnel_plot.pdf", width=7, height=6.5)
mr_funnel_plot(singlesnp_results = res_single)
dev.off()

左右两边snp相对对称的话，说明数据没有异质性。

# 15.留一法敏感性分析
pdf(file="leaveoneout.pdf", width=7, height=6.5)
mr_leaveoneout_plot(leaveoneout_results = mr_leaveoneout(dat))
dev.off()

使用留一法每次去除一个snp然后计算效应值和波动范围。然后还可以把去除某一个snp后的值与所有snp的效应值进行比较，如果波动情况不大的话就说明去除单个snp对孟德尔随机化结果没有产生过多的影响。

参考资料：

TwoSampleMR： https://github.com/MRCIEU/TwoSampleMR

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台。更多相关内容可关注公众号：生信方舟 。

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

数据挖掘

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

数据挖掘

登录后参与评论

0 条评论

热度