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
-
生成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 -
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_esiwrite.table(eqtl_esi,"eqtl_1.esi",quote = F,col.names = F,row.names = F)
更新esi
smr --beqtl-summary eqtl --update-esi eqtl_1.esi
