当孟德尔遇上GWAS…

Hi，大家好，我是晨曦

今天这期推文是我们孟德尔随机化分析系列推文的新起点，前两期推文我们已经大致了解了孟德尔随机化分析（MR分析）的总体流程，尽管整体的分析流程比较固化，但是其中的一些细节还有难点依旧值得我们讨论和学习，如果在进行MR分析的时候遇到了以下问题，那么请关注随后几期的MR分析教程，相信会让各位小伙伴对于MR分析有更深入的理解：

1.从哪里获取数据？

2.使用什么R包进行分析？

3.工具变量的选择有什么技巧？

4.混杂因素如何校正，以及在文章书写的时候究竟需要如何描述？

..........

好，那么本期推文是MR分析系列教程的第一期，我们终点关注如何获取数据，因为在MR分析中数据是核心，如果我们没有数据那么其实一切都是纸上谈兵|镜花水月而已

那么，我们开始吧~

一、MR分析数据库介绍

首先推荐给大家的第一个GWAS数据库就是——IEU OpenGWAS project (mrcieu.ac.uk)

优点：可以和TwoSampleMR包完美配合，十分方便

缺点：数据部分时间有点old，而且大部分我们感兴趣的数据其实不能通过R包进行完美下载

那么，第二个GWAS数据库是——UK Biobank — Neale lab

优点：数据全面、整理有条理

缺点：UKbiobank是英国数据库，但是貌似完整版是需要花钱注册的（UKB数据使用费用分为三档，分别是3000英镑/6000英镑/9000英镑），我们这里使用的其实是2018年有UKbiobank分享出来的数据，这个数据是免费且不需要注册的，所以数据时间上可能会有点滞后，但是也还好，因为不像是测序数据，一般来说大批量的GWAS数据普通课题组是无法承担的

随后，第三个GWAS数据库是——FinnGen-tutkimushanke vie suomalaiset löytöretkelle genomitietoon | FinnGen

优点：数据较新、包括新型冠状病毒的GWAS数据、而且页面布局超级好看！！

当然还有很多的GWAS数据库，请注意，MR分析其实依靠的就是GWAS数据，但是GWAS数据并不是只有MR分析这一种分析思路，所以当我们掌握了MR分析之后，也可以去开展其它的GWAS分析流程

而且网络还有很多总结好的GWAS数据库，这里晨曦只是展示了自己比较常用的，更多的数据库，各位小伙伴可以在网络上进行检索就可以找到更多总结好的教程

二、数据读取

那么，下面我们通过IEU数据库来展示两种数据读取的方式，使用IEU数据库来进行演示，一方面这个数据库有配套的R包可以做到完美配合，一方面就是这个数据库上也有不能通过配套R包进行完美配合，但是可以下载的GWAS数据

其实简单来说就是在线方式和本地方式两种数据获取方式的演示

本地方式

假设：BMI和急性心肌梗死是否具有因果关系，即BMI升高会不会导致急性心肌梗死的发病率增加

这个主题会一直持续到本次系列推文结束，我们这期主要涉及到获取数据的内容

第一步：获取暴露因素的GWAS数据，打开IEU数据库搜索相关GWAS数据

直接在IEU的搜索栏中搜索相关疾病即可，这里有有一个细节，就是我们搜索的关键词尽可能不要太狭窄，其实就是和文献检索一个道理，先保证查全，再靠人工实现查准

然后我们就可以看到上面这个界面，这个时候需要注意GWAS ID，我们可以很清楚的知道，上面这些GWAS数据都是来自于UKbiobank的，所以这些数据我们无法通过配套的R包直接获得，然后我们随便点击一个GWAS ID

然后重点看我标志箭头的地方，VCF文件就是我们需要下载的GWAS全部数据文件，然后这个数据的时间是2018年，这里晨曦想要提出一个疑问：为什么是2018年呢？

回答：因为目前网络上，免费，开源，贡献的UKbiobank数据都是截至在2018年的，后续更新的数据需要注册且缴纳一定费用

然后这个数据是欧洲人的数据，其中心肌梗死的患者有1294人，正常健康人有461716人，这就是为什么GWAS数据我们一般只能靠挖掘，就是因为数据量是很庞大的

这里晨曦建议，如果我们暴露因素选择是欧洲人，那么我们的终点结局也需要选择欧洲人，尽可能做到人种的统一，同时我们的疾病组尽量保证在大于1000人以上，如果太少，除非必须，否则晨曦就不建议使用这个GWAS数据了

三、数据信息

然后这里我们点击Download VCF就可以下载这个数据集，然后我们通过下面代码就可以在R中打开

#准备工作library(TwoSampleMR)library(tidyverse)library(vcfR)#读取VCF文件data <- vcfR::read.vcfR("ukb-b-453.vcf.gz")data

然后我们可以看到这个VCF对象里面有三部分数据，然后我们这里简单查看一下每一个部分的数据

meta <- data.frame(data@meta)meta

然后我们可以知道meta信息里存放的是注释信息，也就是每一个变量是什么意思

在前面的推文中我们提到过，我们想要进行MR分析需要必须的信息包括：

1.SNP ID ：通常是rs开头，但是如果遇到不是的时候，我们可以通过查询网络看看可不可以互换

2.Effect allele（alternative allele）：效应位点

3.Other allele（reference allele）：其它位点

4.Beta（OR）：如果是连续型变量就是Beta，如果是二分类变量就是OR，代表突变究竟有益还是有害

5.Standard error：置信区间

6.Pvalue ：SNP位点是否具有统计学意义

那么很显然这个meta信息就是帮助我们理解变量究竟是什么含义的，因为有可能变量代表某个含义但是写法不一样，所以需要单独注释出来，这就是meta信息存在的意义

fix <-data.frame(data@fix)fix

很显然，fix信息就是SNP矩阵，但是有些信息却是我们不需要的，我们还需要查看剩下的信息

REF：reference allele（Other allele）

ALT：alternative allele（Effect allele）

gt <- data.frame(data@gt)gt

gt信息包含的就是一些统计学的信息，通过查询meta信息我们可以知道每一个信息的含义

ES：Effect size（其实就是beta值）

SE：置信区间（标准误）

LP：-log10 p-value

AF：EAF（效应位点突变频率）

四、数据处理与分析

然后既然我们已经明确了上述信息都是代表什么，那么我们就可以把数据进行合并和删除，整理成MR分析需要的样子

肯定有很多小伙伴不明白这些信息具体是什么意思，这个因为涉及到遗传学的一些知识，而且如果一定要问晨曦这些变量的具体含义，晨曦也不是十分的理解，但是我们只需要知道我们需要这些信息来进行后续的分析，如果真的没有到时候在进行相关的网络搜索即可

#读取数据data <- vcfR::read.vcfR("ukb-b-453.vcf.gz")meta <- data.frame(data@meta)fix <- data.frame(data@fix)gt <- data.frame(data@gt)#整理数据fix <- data.frame(data@fix[,1:5])fix <- fix %>% dplyr::select(ID,ALT,REF,everything())gt <- data.frame(data@gt[,2])beta <- as.numeric(unlist(strsplit(as.character(gt$data.gt...2.),split = ":"))[seq(1,nrow(gt)*5,5)])se <- as.numeric(unlist(strsplit(as.character(gt$data.gt...2.),split = ":"))[seq(2,nrow(gt)*5,5)])p <- as.numeric(unlist(strsplit(as.character(gt$data.gt...2.),split = ":"))[seq(3,nrow(gt)*5,5)])MR_data <- data.frame(beta = beta, se = se, adjpvalue = p)MR_data$pvalue <- 10^(MR_data$adjpvalue)MR <- cbind(fix,MR_data)rownames(MR) <- NULLcolnames(MR)[1] <- "SNP"colnames(MR)[c(2,3)] <- c("Effect_allele","Other_allele")MR#write.csv(MR, file="exposure.csv")

然后这样我们就得到了暴露因素和工具变量之间的GWAS数据，然后我们需要选择工具变量和暴露因素具有强相关性，那么原始P值需要小于5×10-8且LD<0.001，最好在计算一个F大于10（但是有的文献并没有计算这个，而是单纯拿前面两条进行筛选），关于工具变量的选择条件（MR分析的三大假设）会专门出一期推文来进行讲解

那么到这里我们通过下载数据并整理成R包需要的方式就结束了，然后我们需要使用TwoSampleMR包进行分析

myocardial_exposure_clump <-read_exposure_data(filename = "exposure.csv",sep = ",",snp_col = "SNP",beta_col = "beta",se_col = "se",effect_allele_col = "Effect_allele",other_allele_col = "Other_allele",clump = TRUE)