Make a BESD file from QTLtools output

Make a BESD file from QTLtools output

# 根据基因名与变异rsid编号去重   
library(tidyverse)
library(readxl)
qtltoolsnomi <- read_table("qtltoolsnomi.txt",col_names = F)
qtltoolsnomi %>% distinct(X1,X8, .keep_all = TRUE) -> qtltool
qtltool %>% group_by(X2) %>% 
  summarise(n=n()) %>% 
  as.data.frame()
# 去除X,Y染色体
qtltool %>% filter(!(X2 %in% c("chrX","chrY"))) -> qtltool
qtltool %>% filter(!(X8==".")) -> qtltool
write.table(qtltool, "qtltool.txt",quote = FALSE,row.names = FALSE, col.names = F)

生成 BESD file

# 工作目录
$ pwd
/HP_NL/PUBLIC/WORK/guanjiawei/analysis/TWAS_papper/eQTL_result
$/HP_NL/PUBLIC/WORK/guanjiawei/bio_software/soft/smr/smr --eqtl-summary qtltool.txt --qtltools-nominal-format --make-besd --out eqtl_ali_225

Update MAF file to a BESD file

1. 生成MAF文件

   
   awk '{print $2,$3,$4, $5}' ./snp_info.frq > ./eqtl.freq
   # 删除第一行列名
   $ sed -i '1d' eqtl.freq
   
   # 去除包含NA的行
   grep -v 'NA' eqtl > a.txt
   mv a.txt  eqtl.freq
   
   $ head eqtl.freq 
   rs3107975 C T 0.1477
   rs76735897 G A 0.175
   rs77573425 C G 0.1727
   

2. update the frequencies of the effect alleles，实质更新到.esi文档

  library(tidyverse)
   library(dplyr)
   eqtl_esi <- read_table("eqt.esi",col_names = F)
   names(eqtl_esi)[2] <- "rsid"
   eqtl_freq <- read_table("eqtl.freq",col_names = F)
   names(eqtl_freq) <- c("rsid","A1","A2","freq")
   eqtl_esi %>% 
     left_join(eqtl_freq, by = "rsid") %>% 
     dplyr::select(X1,rsid,X3,X4,A1,A2,freq) -> eqtl_esi
   write.table(eqtl_esi,"eqtl_1.esi",quote = F,col.names = F,row.names = F)

更新esi

smr --beqtl-summary eqtl --update-esi eqtl_1.esi