CAUSE 使用记录

onlooker_zsq

已于 2023-06-01 17:23:41 修改

阅读量1.1k

点赞数 2

文章标签： r语言

于 2022-09-30 14:10:16 首次发布

本文链接：https://blog.csdn.net/onlookerq/article/details/127122587

版权

library(readr)
library(dplyr)
library(cause)
library(ieugwasr)
X1=read.table('4_IBS.txt',sep = '\t',header = T)
X2=read.table('PTSD_2019.txt',sep = '\t',header = T)
X=gwas_merge(X1,X2,snp_name_cols = c("SNP","SNP"), 
                  beta_hat_cols = c("b","BETA"), 
                  se_cols = c("se","SE"), 
                  A1_cols = c("A1","A1"), 
                  A2_cols = c("A2","A2"),
                  pval_cols = c('p','P')
             )
set.seed(100)
varlist <- with(X, sample(snp, size=1000000, replace=FALSE))
params <- est_cause_params(X, varlist)
params$rho
r2_thresh = 0.01
pval_thresh = 1e-3
X_clump <- X %>%
  rename(rsid = snp,
         pval = p1) %>%
  ieugwasr::ld_clump(dat = .,
                     clump_r2 = r2_thresh,
                     clump_p = pval_thresh,
                     plink_bin = genetics.binaRies::get_plink_binary(), 
                     pop = "EUR",
                     bfile = './EUR/EUR')
#windows系统下plink不匹配问题，直接去plink官网下载相应版本放入对应的文件夹中，若出现无法发现目标文件的错误，进入目标文件的文件夹，将文件夹中所有内容均删除之后再尝试
top_vars <- X_clump$rsid
res <- cause(X=X, variants = top_vars, param_ests = params)
res$elpd
summary(res, ci_size=0.95)$tab
summary(res)$p
summary(res)$tab
plot(res)
elpd_table <- recompute_elpd_table(res)
pnorm(res$elpd$z[3])

ld_clump使用本地EUR文档，ieugwasr::ld_clump功能

详细教程见Example Analysis with CAUSE: LDL -> CAD