TwoSampleMR:local clump(MR-Base exceeded 300 seconds) 包括Windows和Linux R解决办法

置顶 Bioinfoer

已于 2024-06-12 19:10:40 修改

阅读量1.9w

点赞数 29

文章标签： r语言

于 2022-08-10 20:04:20 首次发布

本文链接：https://blog.csdn.net/xiaozheng1213/article/details/126269969

版权

首先是Windows

一个做孟德尔随机化的过程遇到的报错：

bmi_exp_dat <- clump_data(bmi_exp_dat,clump_r2=0.01,pop = "EUR")

Please look at vignettes for options on running this locally if you need to run many instances of this command.
Clumping C5nTuK， 5340156 variants, using EUR population reference
Error in api_query("ld/clump", query = list(rsid = dat[["rsid"]], pval = dat[["pval"]], :
The query to MR-Base exceeded 300 seconds and timed out. Please simplify the query

以下是作者给出的解决办法：

Extracting without clumping
Performing clumping on each chromosome separately

Alternatively, you could try downloading the VCF from IEU OpenGWAS project and using the gwasvcf to extract based on p-value, and use the ieugwasr package to do clumping locally. We're trying to get to a point where it's easy to do heavier computation on these data locally.

LD reference files listed in the gwasvcf page

看到作者建议可以用ieugwasr这个包做本地clump。

下面尝试做本地clump

首先安装ieugwasr包

devtools::install_github("mrcieu/ieugwasr")

然后我们可以看到

ld_clump_local()

Wrapper for clump function using local plink binary and ld reference dataset

ld_clump_local()这个函数是用来做本地clump的，用法和参数如下：

ld_clump_local(dat, clump_kb, clump_r2, clump_p, bfile, plink_bin)

Arguments

dat	Dataframe. Must have a variant name column ("variant") and pval column called "pval". If id is present then clumping will be done per unique id.
clump_kb	Clumping kb window. Default is very strict, 10000
clump_r2	Clumping r2 threshold. Default is very strict, 0.001
clump_p	Clumping sig level for index variants. Default = 1 (i.e. no threshold)
bfile	If this is provided then will use the API. Default = NULL
plink_bin	Specify path to plink binary. Default = NULL. See https://github.com/explodecomputer/plinkbinr for convenient access to plink binaries

需要plink包，安装一下plink包，获得可执行二进制plink地址

devtools::install_github("explodecomputer/plinkbinr")
library(plinkbinr)
get_plink_exe()
#[1] "D:/R-4.1.1/library/plinkbinr/bin/plink_Windows.exe"

然后我们下载bfile，我的GWAS是欧洲人种所以用EUR*

wget http://fileserve.mrcieu.ac.uk/ld/1kg.v3.tgz

#ld_clump_local()报错，作者是把两个函数合并了，在ld_clump函数里包含bfile参数的话，会自动调用ld_clump_local

所以改用ld_clump

ld_clump( dplyr::tibble(rsid=dat$rsid, pval=dat$pval, id=dat$trait_id), plink_bin = genetics.binaRies::get_plink_binary(), bfile = "/path/to/reference/EUR" )"

b <- ld_clump(
    dplyr::tibble(rsid=a$rsid, pval=a$p, id=a$id),
    #get_plink_exe()
    plink_bin = "D:/R-4.1.1/library/plinkbinr/bin/plink_Windows.exe",
    #欧洲人群参考基因组位置
    bfile = "D:/EUR_ref/EUR"
)

注意a的列名，必须要有:（with the following columns:）

rsid
pval
trait_id

然后就完成了，可以看到跟在线方法去除的LD SNP个数是一样的

PLINK v1.90b6.10 64-bit (17 Jun 2019)          www.cog-genomics.org/plink/1.9/
(C) 2005-2019 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to C:\Users\Lenovo\AppData\Local\Temp\RtmpS8aq1j\file2e646a65362c.log.
Options in effect:
  --bfile D:/EUR_ref/EUR
  --clump C:\Users\Lenovo\AppData\Local\Temp\RtmpS8aq1j\file2e646a65362c
  --clump-kb 10000
  --clump-p1 0.99
  --clump-r2 0.001
  --out C:\Users\Lenovo\AppData\Local\Temp\RtmpS8aq1j\file2e646a65362c

32549 MB RAM detected; reserving 16274 MB for main workspace.
8550156 variants loaded from .bim file.
503 people (0 males, 0 females, 503 ambiguous) loaded from .fam.
Ambiguous sex IDs written to
C:\Users\Lenovo\AppData\Local\Temp\RtmpS8aq1j\file2e646a65362c.nosex .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 503 founders and 0 nonfounders present.
Calculating allele frequencies... done.
8550156 variants and 503 people pass filters and QC.
Note: No phenotypes present.
--clump: 78 clumps formed from 2011 top variants.
Results written to
C:\Users\Lenovo\AppData\Local\Temp\RtmpS8aq1j\file2e646a65362c.clumped .
Warning: 'rs9930333' is missing from the main dataset, and is a top variant.
Warning: 'rs8083289' is missing from the main dataset, and is a top variant.
Warning: 'rs2683992' is missing from the main dataset, and is a top variant.
27 more top variant IDs missing; see log file.
Removing 1963 of 2041 variants due to LD with other variants or absence from LD reference panel

Linux同理

expo_dat <- expo_dat[which(expo_dat$pval.exposure<1e-5),]

b <- ld_clump(
    dplyr::tibble(rsid=expo_dat$SNP, pval=expo_dat$pval.exposure, id=expo_dat$id.exposure),
    #plink位置
    plink_bin = "/GM_GWAS_LD_clumped_snps/plink",
    bfile = "/GM_GWAS_LD_clumped_snps/EUR_ref/EUR",
    clump_kb = 1000,clump_r2 = 0.1
)
expo_dat <- expo_dat[which(expo_dat$SNP %in% b$rsid),]

Clumping A9o7Sb, 50 variants, using EUR population reference
PLINK v1.90b6.21 64-bit (19 Oct 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to /tmp/RtmpdrrA6m/file4bc2f27623281.log.
Options in effect:
  --bfile /share2/pub/zhenggw/zhenggw/GM_GWAS_LD_clumped_snps/EUR_ref/EUR
  --clump /tmp/RtmpdrrA6m/file4bc2f27623281
  --clump-kb 1000
  --clump-p1 0.99
  --clump-r2 0.1
  --out /tmp/RtmpdrrA6m/file4bc2f27623281

515328 MB RAM detected; reserving 257664 MB for main workspace.
Allocated 61144 MB successfully, after larger attempt(s) failed.
8550156 variants loaded from .bim file.
503 people (0 males, 0 females, 503 ambiguous) loaded from .fam.
Ambiguous sex IDs written to /tmp/RtmpdrrA6m/file4bc2f27623281.nosex .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 503 founders and 0 nonfounders present.
Calculating allele frequencies... done.
8550156 variants and 503 people pass filters and QC.
Note: No phenotypes present.
Warning: 'rs55654746' is missing from the main dataset, and is a top variant.
--clump: 8 clumps formed from 49 top variants.
Results written to /tmp/RtmpdrrA6m/file4bc2f27623281.clumped .
Removing 42 of 50 variants due to LD with other variants or absence from LD reference panel

更新了更新了

新增加In file(file, "rt") :
cannot open file '/tmp/Rtmp8ur92k/xxx': No such file or directory 解决办法

好多人出现了这个问题，我今天做LD matrix也遇到了，来解决一下，报错信息是找不到文件，我考虑是tmp文件在程序运行过程中被自动删除了。尝试手动添加tmp文件名

ld_matrix_local <- function(variants, bfile, plink_bin, with_alleles=TRUE)
{
	# Make textfile
	shell <- ifelse(Sys.info()['sysname'] == "Windows", "cmd", "sh")
	fn <- tempfile()


ld_matrix_local <- function(variants, bfile, plink_bin, with_alleles=TRUE)
{
    # Make textfile
    shell <- ifelse(Sys.info()['sysname'] == "Windows", "cmd", "sh")
    fn <- "~/plinkLD/tmpfile/tmp"

发现还是同样error

但是提示里面发现Using up to 39 threads (change this with --threads)，考虑可能是线程冲突，导致文件没有产生，所以在源代码加入线程变量

    fun2 <- paste0(
        shQuote(plink_bin, type=shell),
        " --bfile ", shQuote(bfile, type=shell),
        " --extract ", shQuote(fn, type=shell), 
        " --r square ", 
        " --keep-allele-order ",


        " --threads 1 ",


        " --out ", shQuote(fn, type=shell)
    )

515328 MB RAM detected; reserving 257664 MB for main workspace.
Allocated 61144 MB successfully, after larger attempt(s) failed.
8550156 variants loaded from .bim file.
503 people (0 males, 0 females, 503 ambiguous) loaded from .fam.
Ambiguous sex IDs written to /tmp/Rtmp8ur92k/file557ea8a4015e.nosex .
--extract: 3839 variants remaining.
Using 1 thread.
Before main variant filters, 503 founders and 0 nonfounders present.
Calculating allele frequencies... done.
3839 variants and 503 people pass filters and QC.
Note: No phenotypes present.
--r square to /tmp/Rtmp8ur92k/file557ea8a4015e.ld ... done.

顺利解决！！！

所以原函数加一行代码就好 --threads 1 ，这个线程1好像速度也蛮快的，不懂有什么影响