批量对显著性SNP进行注释：bedtools

大家好，我是邓飞。

之前介绍过使用bedtools进行基因注释的方法：

snpEFF和bedtools基因注释

使用bedtools进行gwas基因注释

虽然，饭要一口一口的吃，注释要一个一个的进行，否则走的太快……。然鹅，不批量彰显不出数据分析师的骄傲。

本来，同一个性状的显著性位点，一起进行注释，这，不算什么。但是多个性状的显著性位点，放在一起，批量注释，而且还可以区分出不同性状的结果，这，就需要一点技术了，这，也是今天博客的目的！

1，处理显著性的SNP

将位点按照：性状，染色体，物理位置进行整理

2，整理gff文件

默认的官方gff文件包括九列，。注意类型中匹配gene

awk -F "\t" 'BEGIN {OFS="\t"}$3=="gene" {print $1,$4,$5,$9}' genes.gff > Chr_position_gene

有时候染色体编号不对应，需要重新编码对应起来。

格式为gff3就行，注意分隔符需要是tab，不能是空格！

3，将不同性状的显著SNP保存为txt文件

nn = table(snp$Trait) %>% names
fwrite(snp[snp$Trait==nn[1],2:3],paste0(nn[1],".txt"),col.names = F,quote = F,sep = " ")fwrite(snp[snp$Trait==nn[2],2:3],paste0(nn[2],".txt"),col.names = F,quote = F,sep = " ")fwrite(snp[snp$Trait==nn[3],2:3],paste0(nn[3],".txt"),col.names = F,quote = F,sep = " ")fwrite(snp[snp$Trait==nn[4],2:3],paste0(nn[4],".txt"),col.names = F,quote = F,sep = " ")fwrite(snp[snp$Trait==nn[5],2:3],paste0(nn[5],".txt"),col.names = F,quote = F,sep = " ")fwrite(snp[snp$Trait==nn[6],2:3],paste0(nn[6],".txt"),col.names = F,quote = F,sep = " ")fwrite(snp[snp$Trait==nn[7],2:3],paste0(nn[7],".txt"),col.names = F,quote = F,sep = " ")fwrite(snp[snp$Trait==nn[8],2:3],paste0(nn[8],".txt"),col.names = F,quote = F,sep = " ")fwrite(snp[snp$Trait==nn[9],2:3],paste0(nn[9],".txt"),col.names = F,quote = F,sep = " ")fwrite(snp[snp$Trait==nn[10],2:3],paste0(nn[10],".txt"),col.names = F,quote = F,sep = " ")

4，将SNP的信息，增加区间

add_qujian.R FLL.txt 50000

代码主要是自动添加显著位点的上下游，并保存，格式用bedtools可以直接处理。

$ cat add_qujian.R #! /home/kitty/bin/Rscriptargs = commandArgs(T)if(length(args) < 2){  cat("\n请添加参数，运行示例：add_qujian.R b 50000，这里b两列染色体和物理位置的显著snp，50000是50k的区间\n\n")  quit("no")}

file = args[1]nn = as.numeric(args[2])
# file = "sig_snp.txt"# nn = 50000
library(data.table)library(tidyverse)
dd = fread(file,header=F)%>% select(Chr=1,V2=2)head(dd)
dd1 = dd %>% mutate(start = V2 - nn, end = V2 + nn) %>% select(Chr=1,start,end)dd1 = dd1 %>% arrange(Chr,start,end)head(dd1)summary(dd1)dd1[dd1$start<0]$start = 1summary(dd1)fwrite(dd1,"snp_re_qujian_sort.txt",sep = "\t",col.names = F,quote = F)

5，写个脚本批量处理

#!/usr/bin/bash
if [ $# != 3 ]then    echo bash $0 sig_snp.txt 5000 t2.gff3; 注意，sig_snp.txt是染色体和物理位置两列，50000是50k的上下游区间，gff3是注释文件    exit 1fi
add_qujian.R sig_snp.txt $2bedtools merge -i snp_re_qujian_sort.txt >snp_qujian.bed#nn=${3##*/}bedtools intersect -a $3 -b snp_qujian.bed -wa >${1##*/}_sig_snp_plus_gene_result.txt

6，用xargs循环批量处理

echo *txt|xargs -n1|xargs -i bash run_bedtools_from_sig_snp.sh {} tt.gff3

飞哥后记：今天没有什么素材，把自己之前总结的代码搞出来，发现也还不错。希望对大家有所帮助，后面这些内容都会连同数据代码放到我的GWAS cookbook中，下一版的更新中吧。