A package for running MR In batches and in parallel quickly

由于个人课题安排，时间紧张以及盗版猖獗等原因，get_mr暂时停止更新 已失效的网盘链接已重新更新如下：https://pan.baidu.com/s/1yyYA4cVHKnyCce9YSPRA9A?pwd=ipvl 提取码：ipvl 万望大家理解，感谢陪伴

1. 禁止一切倒卖我们代码的行为。我们承诺已开源代码永久性免费开源，公开可得
2. 使用代码前，务必确认代码无误。我们没有利用这个R包发表文献，也没有进行倒卖，或授课。因此，我们无法确保我们的代码每一环，完全的，彻底的，100%的，没有错误。我们的代码本质上是分享经验，我们并没有因为这些代码获得任何利益，也没有利用它来当成文章发表，虽然我们会尽力，但我们无法保证彻底的完美。我们对所有因代码可能存在的错误导致的纠纷不负有任何责任。
3. 我们不会参与任何商业行为和组织，我们没有与任何组织进行合作。任何倒卖行为都是侵权行为，非倒卖的利用行为属于第三方的行为，与我们无关，由第三方对本身内容负责。因代码错误导致的纠纷参考第2条声明。
4. 我们分批开源不是利益因素导致的，我们没有将更新的版本的Get_MR进行售卖！我们本意是为了保护内部小伙伴的努力成果，优先使用我们最新开发的前沿功能，万望大家理解！

1. 修复cyclemr函数中clean函数bug
2. 补充说明，get_outcome函数修改了TwoSampleMR包原函数的默认值，将maf阈值调整为0.4，并不获取proxy（默认值为0.3与获取proxy）。由于这个函数主要是为了方便批量预实验（添加了错误防停止运行机制），因此建议正式做数据采用官方函数。

欢迎来到向量化与并行化的世界 本次更新就一个重大功能，就是并行化运行mr分析，以最优的效率批量跑大量的数据。家用电脑较新服务器基本可以实现2小时批量运行10000个因素，如果你有高性能服务器，恭喜你，30分钟以内就能跑完。

由于本次主要是思路与方法的分享，所以函数的帮助文档写的不多，主要还是看示例代码即可，应该还是很容易上手的。

公众号回复：“示例”即得示例代码 (不会还没关注我们公众号吧#doge，GetScience, 等你来！)

这个函数是用于执行循mr 分析的功能函数，可以在R语言中使用。该函数可以将数据分配到多个计算节点中运行，提高MR分析的效率。

直接在R语言中调用这个函数，如下所示：

# 调用cyclemr函数
mr_results <- cyclemr(dat = data, cl_num = 4, type = "list")

• dat: harmonise_data后的数据，可以是数据框或列表类型。默认是list
• cl_num: 批量化线程数。
• type: MR分析数据类型，可以是"list"或"data"，默认为"list"。主要使用情况也是list

注： 如何判断自己的电脑能开启的最大线程数？

在任务管理器可看到

• 内核： 计算机核心数

• 逻辑处理器： 计算机总线程数

  比如说很多处理器宣传是8核16线程，这个8就指的是内核，这个16指的是逻辑处理器。

  本质上，16只是将每个核心一分为2，但是他们能干的活是一样的。所以一般设置为内核数即可满载CPU

  当然这不能一概而论，因为每个CPU和厂家调度不一样，如果你发现使用内核数不能让CPU跑到100%，则尝试用逻辑处理器数

• cyclemr函数返回一个包含MR分析结果的数据框。

下面是使用这个函数的一个示例：

mr_results <- cyclemr(dat = data, cl_num = 16, type = "list")

在设置无误情况下，这是我手头有的所有电脑测试出的运行时间：

运行10000个数据。（ieu批量数据前10000个）除了服务器外，其他均使用Windows系统

欢迎各位补充自己手头的机器的运行时间数据，尤其M系列的苹果处理器数据

根据CHR和POS，从ensemble官网中获取rsID。

get_rsid(chr, pos, version = 'hg38')

• chr：染色体号。
• pos：基因位置。
• version：表示使用的基因组版本，默认为最新的版本（'hg38'）。也可选择hg19.
• 注： GRCh37=hg19，GRCh38=hg38

该函数基于生物信息学数据库Ensembl SNP Mart来查询给定位置的相关信息。如果未指定基因组版本，则默认使用最新的版本(hg38)。该函数会根据指定的基因组版本选择正确的URL。如果您想查询其他版本的数据，可以将version参数设置为相应版本的字符串。

参考：How to find rsID with biomaRt in R (bioconductor.org)

ds4 <- data.frame(CHR = c("8", "8", "8", "8", "8"),POS = c('101592213', '106973048', '108690829', '102569817', '108580746'))
res<-get_rsid(chr=ds4$CHR, pos=ds4$POS, version = 'hg38')

如果出现：

Error in curl::curl_fetch_memory(url, handle = handle) : 
  Timeout was reached: [grch37.ensembl.org:80] Operation timed out after 300013 milliseconds with 7909 bytes received

这并不是代码问题，而是网络超时了，ensemble的API经常拥堵，多试几次即可。当然也有可能请求的数据量太大，也可能会出现这个问题。

这两个函数是用于进行双样本MR（Two-sample Mendelian Randomization）分析的数据处理和提取过程的功能函数。

• get_exposure函数用于从给定的遗传仪器ID中提取出暴露（exposure）数据。
• get_outcome函数用于从给定的遗传仪器ID和暴露数据中提取出结果（outcome）数据。

使用这两个函数前，需要先安装并加载TwoSampleMR包。

library(TwoSampleMR)

然后可以直接在R语言中调用这两个函数，如下所示：

# 调用get_exposure函数
exposure_data <- get_exposure(id = "ieu-a-1", pval = 5e-8, r2 = 0.001, kb = 10000)

# 调用get_outcome函数
outcome_data <- get_outcome(id = "ieu-a-1", expo = exposure_data)

• get_exposure函数的参数说明：
  • id: 遗传仪器ID。
  • pval: 提取暴露数据的P值阈值，默认为5e-08。
  • r2: 遗传仪器间的LD（linkage disequilibrium）值的阈值，默认为0.001。
  • kb: 遗传仪器的范围（以kb为单位），默认为10000。
• get_outcome函数的参数说明：
  • id: 遗传仪器ID。
  • expo: 暴露数据，TwoSampleMR包格式

获取OPEN GWAS数据库所有可用的ID。可以指定获取某个前缀的ID

使用这个函数前，需要先安装并加载TwoSampleMR包。然后可以直接在R语言中调用这个函数，如下所示：

# 调用get_ao函数
ao <- get_ao()## 不限定来源则返回所有id
ao <- get_ao(a = "finn")## 这样就会返回finn的所有可用id

• a: （可选）数据来源。

来自OPEN GWAS Browse the IEU OpenGWAS project (mrcieu.ac.uk) 5.1获取

（主要用于向量化）用于清洗outcome。将exposure中（list形式，也就是批量化的形式存在的exposure）不存在的SNP从outcome中剔除，大幅精简outcome，并大幅提升harmonise_data的速度。

实测清洗与不清洗outcome对比，速度相差一百倍以上。

直接在R语言中调用这个函数，如下所示：

# 调用clean_outcome_from_exposure函数
cleaned_outcome <- clean_outcome_from_exposure(expo = exposure_data, outcome = outcome_data)

• expo: 暴露数据，格式为list，TwoSampleMR包格式。比如get_exposure批量化获取下来的数据
• outcome: 结果数据，TwoSampleMR包格式。

这个函数是用于清洗遗传关联数据集，使其符合特定的数据集要求的功能函数。

直接在R语言中调用这个函数，如下所示：

# 调用clean_GWAS函数
cleaned_GWAS_list <- clean_GWAS(list = GWAS_list, clean = c("bbj", "eqtl"))

• list:一个list。里面包含每个暴露的data.frame。 具体参考批量运行get_exposure后的结果
• clean: 需要清洗的数据集名称，一个字符型向量类型。具体参考get_ao的附注。

• clean_GWAS函数返回一个清洗后的遗传关联数据集列表，符合特定数据集要求。

4.22：（重要） 修复PRESSO的bug，请使用PRESSO时务必更新代码！！

4.19：注意TwoSampleMR包有一列是mr_keep。请在运行我们get_mr包的所有分析函数前，如果需要使用TwoSampleMR包格式的数据，请将mr_keep状态为false的删掉（如果是没有harmonise之前，列名为mr_keep.exposure/mr_keep.outcome）。因为带有未通过质控的SNP，可能会带来未知的错误！可以参考以下代码

  
  dat<-harmonise_data(exposure,outcome)  
  
  ## 在运行完harmonise_data后可能会产生mr_keep状态为false的，不能用于MR分析的SNP，我们把他删掉：
  
  dat<-subset(dat,mr_keep==T)
  
  
  ## 然后再继续使用各种功能

github：HaobinZhou/Get_MR: A package for running MR In batches and in parallel quickly (github.com)

如果觉得好用，可以点一下github项目上的小星星吗，这是我们继续开源的最大动力，谢谢！

R包以R脚本的形式提供，打开R包，全选运行，即得到所有function

1. 进入githubHaobinZhou/Get_MR: A package for running MR In batches and in parallel quickly (github.com)，下载代码zip

2. source("./Get_MR1.0.r") ## 填文件所在地址
   
   ## 或者直接打开R文件，全选代码运行也可以！

1. 本地clump，1000G处理好的MAF文件，MRlap依赖文件如何获取： GetScience公众号可免费获取已处理好文件，回复"依赖文件"即得链接。源文件请查看本文相应function介绍处

2. 输入clump文件路径后总是报错：

   #尤其注意这个文件名的书写，因为他们是二进制文件，不需要写后缀！只需要选取对应的人种即可，比如欧洲人：
   LD_file="S:/GWAS数据/本地LD依赖文件/EUR"
   
   ## 这个问题我回答好多遍啦！

3. 第一次使用如何安装关联R包：Get_MR/1.0 at main · HaobinZhou/Get_MR (github.com)

   1. 如果不需要使用MungeSumstats包（相关函数包括：format_Mun，get_chr_pos，format_getmr中source="ukb_nosnp") ，则只需要运行Get_MR1.0dependence.R
   2. 如果需要使用MungeSumstats包，则还需运行Install_Reference_Genome.r 这个包括了hg19和hg38的基因组参考文件，总大小达到了5G！如果直接安装失败，在GetScience公众号回复"基因组参考"可得下载链接，并本地安装（推荐）

4. Bug反馈：代码仅由两人编写，难免出现错误。欢迎提交bug到GetScience公众号后台！

5. 感谢所有Get_MR使用的R包作者，是因为他们我们才得以轻松实现这么多复杂的功能。他们都是开源的，因此我们承诺Get_MR将永久免费开源。这意味着使用者可以随意地修改，分发代码，但前提是遵守：

   1.本代码不得用于任何商业或盈利目的

   2.未经代码作者的同意，本代码不得用于任何形式的销售或商业交易

   3.本代码可以在非商业性的科研、学术研究和个人使用的情况下免费使用

   4.在使用本代码并重新打包并向公众发放时，请引用我们的公众号原文

用于计算两个数据框中SNP之间的遗传相关性（rg）。

LDSC_rg(expo, outcome, an, sample_prev = NA, population_prev = NA,
        ld, wld, chr_filter = c(1:22), n_blocks = 200)

• expo: 一个数据框，其中包含一个遗传暴露指标的多个SNP和它们与结果变量的rg。
• outcome: 一个数据框，其中包含一个结果变量的多个SNP和它们与遗传暴露指标的rg。
• an: 它是一个字符串,目前还没有作用（因为我们提供的依赖文件只有eur的，其他人种还没更新）
• sample_prev: 遗传暴露指标的样本流行病学先验患病率。默认为 NA。
• population_prev: 遗传暴露指标的人群流行病学先验患病率。默认为 NA。
• ld: 本地LD依赖文件
• wld: 本地weighted LD 依赖文件
• chr_filter: 一个整数向量，用于指定要使用的染色体。默认为包含1-22的整数向量。
• n_blocks: 用于计算加权LD矩阵的块数。默认为200。

一个具有以下元素的列表：

• rg: 两个数据框中SNP之间的遗传相关性（rg）。
• pval: rg 的双侧P值。
• N_snps: 参与计算rg的SNP数量。

具体用法参照：mr_lap和LDSC_rg示例.r 可通过公众号GetScience回复示例获取文件

mrlap是一种矫正样本重叠后的双样本MR方法。可用于怀疑有样本重叠的数据中。

R包官网：n-mounier/MRlap: R package to perform two-sample Mendelian Randomisation (MR) analyses using (potentially) overlapping samples (github.com)

mr_lap(expo, outcome, ld, hm3, pval, r2, kb, MR_reverse = 1e-03, save_logfiles = F)

• expo: 数据框，为TwoSampleMR包格式的数据
• outcome: 数据框，为TwoSampleMR包格式的数据
• ld: 数据框，本地LD文件路径
• hm3: 数据框，本地HapMap3文件路径
• pval: 数值，MR 工具变量阈值。
• r2: 数值，clump阈值
• kb: 数值，clump阈值
• MR_reverse: 数值，MR 的方向翻转阈值。
• save_logfiles: 逻辑值，是否保存日志文件。

• res: mrlap 结果。

具体用法参照：mr_lap和LDSC_rg示例.r 可通过公众号GetScience回复示例获取文件

一键式执行cause。可批量化执行多暴露对一结局或一暴露对多结局

## 不并行化运行
cause_getmr(expo, outcome, LD_file, r2 = 0.001, kb = 10000, pval = 1e-05)

## 并行化运行
cl<-makeCluster(2) ## 填你想要的并行化的核数，核数越多，需要的运行内存越大
cause_getmr(expo, outcome, LD_file, r2 = 0.001, kb = 10000, pval = 1e-05,cl=cl)

• expo: TwoSampleMR的暴露格式的数据。
• outcome: TwoSampleMR的暴露格式的数据。
• 注意！expo和outcome，可以是data.frame的形式，也可以是一个list（如list[[1:n]]里都包含数据的data.frame)。但不能outcome和expo同时都是list。当expo或outcome，其中一个为list的情况下，是批量运行一对一的cause。比如我读取了10个暴露和1个结局，将10个暴露lapply读取进来就会是一个list。这时候cause_cyclemr 自动运行每个暴露对结局的cause，也就是批量化执行.

## 比如我这里读取3个暴露文件和1和结局文件
id<-c('a.gz','b.gz','c.gz')
expo<-lapply(id,FUN=fread)
outcome<-fread("outcome.gz")
cl=makeCluster(4)## 内存不够的也可以不并行化运行
res<-cause_getmr(expo, outcome, LD_file, r2 = 0.001, kb = 10000, pval = 1e-05,cl=cl)
stopCluster(cl)
## 这样返回的结果就是3个暴露分别对一个结局的cause结果。

• LD_file: 包含LD信息的PLINK文件名。因为需要大批量地clump，在线clump很容易报错，因此我们采用本地clump。需要本地参考文件。下载地址： http://fileserve.mrcieu.ac.uk/ld/1kg.v3.tgz 。 或关注公众号GetScience直接获取。

#尤其注意这个文件名的书写，因为他们是二进制文件，不需要写后缀！只需要选取对应的人种即可，比如欧洲人：
LD_file="S:/GWAS数据/本地LD依赖文件/EUR"

## 这个问题我回答好多遍啦！

• r2: LD的R平方阈值。默认值为0.001。
• kb: LD的距离阈值（以kb为单位）。默认值为10000。
• pval: 用于LD clumping的p值阈值。默认值为1e-05。
• cl: 并行计算的cluster对象。默认值为NULL。在外部使用cl=makeCluster(n),n为你想并行化的核数。注意核数太多不要爆内存了。

cause_cyclemr函数返回cause结果

RAPS_getmr函数执行基于RAPS的MR并返回结果,并画图

expo<-fread('a.gz')
outcome<-fread('b.gz')
expo<-format_data(...)
outcome<-format_data(...)  ## format_data是TwoSampleMR包的函数，格式化。
expo<-pblapply(expo,pval=1,kb=10000,r2=0.001,LD_file=LD_file,FUN=clean_expo) ## 数据很大，建议本地clump，在线很容易报错
dat<-harmonise(expo,outcome)
res<-RAPS_getmr(dat, dir_figure)

• dat: TwoSampleMR包 harmonise_data后输出的数据
• dir_figure: 保存结果图形的目录。

RAPS_getmr函数返回一个包含基于RAPS的MR结果的数据框。

执行MR-PRESSO

mr_Presso(dat, num = 10000)

• dat: 数据框，包含基因表达和疾病风险关联分析的数据。
• num: 整数，模拟数量。

• mr_presso_res: MR-PRESSO 结果。

dat<-harmonise_data(exposure,outcome) ## TwoSampleMR包的harmonise_data函数输出的结果
mr_presso_res <- mr_Presso(dat, num = 10000)

提取 MR-PRESSO 结果中的主要结果

mr_presso_pval(mr_presso_res)

• mr_presso_res: MR-PRESSO 结果。

• mr_presso_main: MR-PRESSO 主要结果。

mr_presso_main <- mr_presso_pval(mr_presso_res) ##mr_Presso输出的结果

根据 MR-PRESSO 分析结果，将离群值剔除，返回剔除离群值后的dat(我一般称为dat_aj， 也就是 adjusted_data)， 可用于后续的IVW等分析。

mr_presso_snp(mr_presso_res, mr_presso_main, dat, type = "list")

• mr_presso_res: MR-PRESSO 结果。
• mr_presso_main: MR-PRESSO 主要结果。
• dat: 数据框或数据框列表，包含基因表达和疾病风险关联分析的数据。
• type: 字符串，输入数据类型。可选值为 "list" 或 "data"。如果是列表形式的（批量化运行后的结果），就是list，如果是普通数据框就是data

过滤后的数据框或数据框列表。

dat<-harmonise_data(exposure,outcome) ## TwoSampleMR包的harmonise_data函数输出的结果
mr_presso_res <- mr_Presso(dat, num = 10000)
mr_presso_main <- mr_presso_pval(mr_presso_res)
data_aj <- mr_presso_snp(mr_presso_res, mr_presso_main, dat, type = "data")

## 用矫正的data可以用于后续的分析，例如重新计算mr
res_aj<-mr(data_aj)

运用MungeSumstats包标准化GWAS 摘要统计数据（包括hg19和hg38转换）。该函数可以将来自Finngen R8和其他来源的 GWAS 摘要统计数据文件清洗为标准的GWAS文件，并可将基因组位置从 ref_genome 转换到 convert_ref_genome。

format_Mun(file, source = "finn_r8", save_path = NULL, lift = F, ref_genome = "hg38", convert_ref_genome = "hg19")

• file：字符向量或数据框，表示要格式化的 GWAS 摘要统计数据文件或数据框。如果输入的是字符向量，则表示文件的路径。如果输入的是数据框，则表示要格式化的数据框。
• source：字符向量，表示输入文件的来源。默认为 "finn_r8"。
• save_path：字符向量，表示格式化文件要保存的路径。默认为 NULL。
• lift：逻辑值，表示是否将基因组位置从 ref_genome 转换到 convert_ref_genome。默认为 F。
• ref_genome：字符向量，表示 GWAS 摘要统计数据文件使用的参考基因组。默认为 "hg38"。
• convert_ref_genome：字符向量，表示要将基因组位置转换到的参考基因组。默认为 "hg19"。

# 从文件中格式化数据
format_Mun("my_sumstats.txt", save_path = "~/formatted_sumstats", lift = F, ref_genome = "hg38")

# 从数据框中格式化数据
my_sumstats_df <- read.csv("my_sumstats.csv")
format_Mun(my_sumstats_df, save_path = "~/formatted_sumstats", lift = F, ref_genome = "hg38")

#格式化数据并升降版本
format_Mun(my_sumstats_df, save_path = "~/formatted_sumstats", lift = T, ref_genome = "hg38", convert_ref_genome = "hg19") ## 从hg38转为hg19

该函数返回格式化的数据框并将其写入磁盘文件。save_path指定保存的位置

预设的快捷格式化 GWAS 摘要统计数据，这个函数用于将来自多个数据来源的 GWAS 摘要统计数据转换为TwoSampleMR 包所需的格式。

format_getmr(data, type = "exposure", source = "finn_r8")

• data：数据框，表示要格式化的 GWAS 摘要统计数据。
• type：字符向量，表示数据类型，可以是 "exposure" 或 "outcome"。默认为 "exposure"。
• source：字符向量，表示数据来源。默认为 "finn_r8"。目前支持的来源有：
  • "finn_r8": Data download - FinnGen Documentation (gitbook.io)
  • "ukb_nosnp": 尼尔数据库（UKB），因为没有rsid，因此需要匹配（已一键完成）。www.nealelab.is/uk-biobank
  • "Mun"： 来自MungeSumstats包格式化后的数据
  • "covid"： COVID19-hg GWAS meta-analyses round 7 (covid19hg.org)
  • "outcome" ： 已经格式化为TwoSampleMR包的“outcome”格式
  • "exposure"：已经格式化为TwoSampleMR包的“exposure”格式
  • "fast_ukb": fastGWA | Yang Lab (westlake.edu.cn)
  • "bac": 2021年肠菌原文数据 Large-scale association analyses identify host factors influencing human gut microbiome composition - PMC (nih.gov)

my_data <- fread("my_data.gz")
format_getmr(my_data, type = "finn_r8", source = "Mun")

该函数返回格式化的数据框。

这个函数用于格式化 GWAS 摘要统计数据中的表型信息，使其符合命名规范，易于保存为文件（例如批量保存计算R2和F值后的文件）。

主要是为了解决，在Windows系统下，保存文件的名称中不能包含特殊字符，例如:,|。

format_trait(list, short = FALSE, short_num = "40")

• list：列表，表示要格式化的 GWAS 摘要统计数据列表。
• short：逻辑值，表示是否要将表型名称缩短。默认为 FALSE。
• short_num：字符向量，表示缩短表型名称的长度。默认为 "40"。

my_list <- list(data1, data2, data3)
format_trait(my_list, short = TRUE, short_num = "20")

该函数返回格式化后的 GWAS 摘要统计数据列表。

这个函数用于从 VCF 文件中读取摘要统计数据。并保存为压缩文件。默认是.gz为后缀的压缩文件。方便下次读取以及节省空间。

这是因为读取VCF文件将消耗大量电脑资源。我们建议批量读取VCF文件后储存为易于读取的压缩包形式。下次读取方便快捷。因此本函数不会直接返回数据框，而是保存为文件

read_vcf_getmr(file_name, nThread = 8, type = ".gz")

• file_name：字符向量，表示要读取的 VCF 文件名。
• nThread：整数，表示要使用的线程数。默认为 8。
• type：字符向量，表示输出文件类型。默认为 ".gz"。

my_file <- "my_file.vcf"
read_vcf_getmr(my_file, nThread = 4, type = ".gz")

该函数没有返回值，而是将读取的数据写入文件。

这个函数用于从文件中读取 GWAS 摘要统计数据。并返回经过P值筛选的文件。一般用于批量读取大量文件时。比如我要批量读取100个暴露数据，每个数据占用运行内存2G。如果100个，则200G，不是一般电脑可以承受。因此每次读取将直接筛选p值，压缩大小

read_easy(file_name, pval = 5e-08)

• file_name：字符向量，表示要读取的文件名。
• pval：数字，表示筛选摘要统计数据的显著性水平。默认为 5e-08。

my_file <- "my_file.csv"
read_easy(my_file, pval = 1e-06)

该函数返回摘要统计数据的数据框。

从1000G的MAF文件中提取EAF并将其与输入数据匹配。

get_eaf_from_1000G(dat, path, type = "exposure")

• dat：一个数据框，为TwoSampleMR包格式的数据
• path：一个字符串，表示包含1000G MAF文件fileFrequency.frq的目录路径。
• type：一个字符串，表示数据是“exposure”（暴露因素）还是“outcome”（结果），默认为“exposure”。

一个数据框，其中包含输入数据的EAF和类型信息（原始、修正或错误）。

该函数将读取1000G MAF文件，然后将其与输入数据进行匹配。对于每个SNP，该函数将检查输入数据中的效应等位基因与1000G中的效应等位基因是否匹配。

如果不匹配，则将EAF设置为NA并将其类型设置为“error”。对于匹配的SNP，EAF将设置为1000G MAF中的值（如果输入数据的效应等位基因是minor allele），或1-MAF（如果输入数据的效应等位基因是major allele），并将其类型设置为“raw”或“corrected”。

如果type参数设置为“outcome”，则函数将使用输入数据中的结果等位基因来查找EAF。

在处理完所有SNP后，该函数将输出一些有关匹配成功、失败以及NA的信息，以及类型信息的说明。

以下是使用该函数的示例：

dat <- get_eaf_from_1000G(dat, "S:/GWAS数据/本地LD依赖文件", type = "exposure")

# 检查输出
head(dat)

fileFrequency.frq为PLINK1.9输出的，根据1000G数据提取的MAF数据

1000G参考文件下载地址：http://fileserve.mrcieu.ac.uk/ld/1kg.v3.tgz

可自行提取MAF数据。或从GetScience公众号中获取已经提取好的fileFrequency.frq文件

该函数利用MungeSumstats包匹配rsid的染色体位置。

get_chr_pos(dat, type = "exposure")

• dat：一个数据框，TwoSampleMR格式
• type：一个字符串，表示要获取SNP染色体位置和参考序列的SNP类型。可选值为"exposure"或"outcome"。

一个数据框，其中包含输入数据框的信息，以及新列chr.exposure或chr.outcome，表示每个SNP的染色体编号。新列pos.exposure或pos.outcome表示每个SNP在染色体上的位置。

该函数使用format_sumstats和format_data函数从1000G项目中获取SNP的染色体位置和参考序列信息。

以下示例演示如何使用get_chr_pos函数：

# 获取曝露变量SNP的染色体位置和参考序列
exposure_chr_pos <- get_chr_pos(dat, type = "exposure")

# 获取结果变量SNP的染色体位置和参考序列
outcome_chr_pos <- get_chr_pos(dat, type = "outcome")

get_f函数计算样本的F统计量并返回F值大于指定阈值的样本数据。返回的数据包括每个SNP的R2和F值

get_f(dat, F_value = 10)

• dat: TwoSampleMR格式，一定要包含eaf.exposure, beta.exposure, se.exposure, 和samplesize.exposure的数据框。
• F_value: 指定的F统计量的阈值，F值大于该阈值的样本将被返回。默认值为10。

注，本计算公式为

F值： R2值：

公式参考文献：

A Multivariate Genome-Wide Association Analysis of 10 LDL Subfractions, and Their Response to Statin Treatment, in 1868 Caucasians - PMC (nih.gov)

Large-scale association analyses identify host factors influencing human gut microbiome composition - PMC (nih.gov)

mr_dircreate_base函数创建基本的目录结构以保存MR分析结果和图形。

mr_dircreate_base(root_dir, project_name, date = NULL)

• root_dir: 保存结果文件夹的根目录。
• project_name: 项目名称，用于在根目录下创建一个以此命名的子目录。
• date: 日期（可选），用于在子目录名称中添加日期以区分不同日期的结果文件夹。默认为NULL，表示不添加日期。

mr_dircreate_base函数返回一个包含结果文件夹路径的列表。

# 创建结果文件夹
res_dir <- mr_dircreate_base("path/to/root/dir", "project_name", date = "20220327")

# 打印结果文件夹路径
print(res_dir)

• root_dir参数指定保存结果文件夹的根目录。
• project_name参数指定项目名称，用于在根目录下创建一个以此命名的子目录。
• date参数（可选）用于在子目录名称中添加日期以区分不同日期的结果文件夹。默认为NULL，表示不添加日期。
• mr_dircreate_base函数将在根目录下创建一个名为project_name的子目录，并在该子目录下创建4个子目录，分别命名为1.figure、2.table、3.figure of sig res和4.snp with Fval。函数返回一个包含结果文件夹路径的列表。

用于快捷筛选工具变量的函数，可执行P值筛选，EAF值筛选，clump。

## 完整可选参数
clean_expo(expo, pval, low_af = 0.5, high_af = 0.5, clump = TRUE, kb = 10000, r2 = 0.001, LD_file = NULL, af_filter = FALSE)

##不提供LD_file则自动在线clump
clean_expo(expo, pval, clump = TRUE, kb = 10000, r2 = 0.001)
##提供则本地clump
clean_expo(expo, pval, clump = TRUE, kb = 10000, r2 = 0.001，LD_file=LD_file)

• expo: 一个数据框，其中包含遗传暴露指标的SNP名称、beta值、标准误、p值和频率。
• pval: 用于筛选遗传暴露指标的p值阈值。p值小于此阈值的SNP将被保留。
• low_af: 频率过滤的下限值。默认为0.5。如果af_filter为TRUE，则只有遗传暴露指标的频率低于此值或高于high_af时，才会被保留。
• high_af: 频率过滤的上限值。默认为0.5。
• clump: 一个逻辑值，指示是否使用PLINK进行SNP聚类。默认为TRUE。
• kb: 聚类的窗口大小（以kb为单位）。默认为10000。
• r2: LD阈值。默认为0.001。
• LD_file: PLINK二进制文件的路径。如果未提供，则默认在线clump。
• af_filter: 一个逻辑值，指示是否启用频率过滤。默认为FALSE。

用于清理列表中元素行数的R包。常用于批量化质量控制。

clean_list(list, nrow = 10)

• list: 一个列表，其中包含多个元素。
• nrow: 用于筛选元素的行数阈值。如果元素的行数小于此阈值，则该元素将被删除。默认为10。

一个列表，其中仅包含行数大于 nrow 的元素。

# 创建一个包含5个数据框的列表，每个数据框包含1-5行
set.seed(123)
lst <- list(data.frame(a = rnorm(1), b = rnorm(1)),
            data.frame(a = rnorm(2), b = rnorm(2)),
            data.frame(a = rnorm(3), b = rnorm(3)),
            data.frame(a = rnorm(4), b = rnorm(4)),
            data.frame(a = rnorm(5), b = rnorm(5)))

# 运行 clean_list 函数，将阈值设置为2
cleaned_lst <- clean_list(lst, nrow = 3)

# 查看清理后的列表
cleaned_lst

用于从一个数据框中清理具有显著的MR反向因果效应P值的IV。

clean_IV_from_outsig(dat, MR_reverse = 1e-03)

• dat: 一个数据框，其中包含每个IV和其与结果变量之间的MR反向因果效应的P值。
• MR_reverse: 用于筛选IV的MR反向因果效应P值阈值。具有P值小于此阈值的IV将被保留。默认为1e-03。

一个数据框，其中包含P值大于 MR_reverse 值的IV（也就是反向不显著的IVs）。

• 代码作者：广州医科大学 第一临床学院 周浩彬 第二临床学院 谢治鑫

• 帮助文档作者： 周浩彬

• 时间：2023/5/1

• 适配版本: Get_MR2.0

• 开源许可证：GPL3.0

• 公众号： GetScience

• 致谢：感谢广州医科大学 第六临床学院 黄覃耀和 南山学院 林子凯在孟德尔随机化概念，代码思路等提供的重要的建设性建议。