批量对显著性SNP进行注释：bedtools
虽然，饭要一口一口的吃，注释要一个一个的进行，否则走的太快……。然鹅，不批量彰显不出数据分析师的骄傲。
本来，同一个性状的显著性位点，一起进行注释，这，不算什么。但是多个性状的显著性位点，放在一起，批量注释，而且还可以区分出不同性状的结果，这，就需要一点技术了，这，也是今天博客的目的！
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/Rscript
args = 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 = 1
summary(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 1
fi
add_qujian.R sig_snp.txt $2
bedtools 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中，下一版的更新中吧。
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2023-11-09，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
腾讯云开发者
