bedtools | 筛选重合区间 注释bed区间

---

前言

通常我们手里有一个bed区间的范围，但只看这个没有什么概念，所以，要对这个bed区间进行注释，看看都是什么基因。

---

1、对文件格式进行调整

首先我们需要下载一个RefGene的文件，里面包含转录本及外显子等诸多信息。

anno_file = 'gs.anno' # ucsc 下载的refgene文件
bed_file = '44.bed' # 目标bed区间 我的是4列，但只用里前三列
target_gene_file = '44_gene_list.txt' # 想要保留的基因名称
out_file = 'gs.anno.tx' # 结果文件  4 列

with open(target_gene_file) as f_gene:
    target_gene_list = f_gene.read().strip().split('\n')
with open(anno_file) as f_anno:
    anno_list = [i.split('\t') for i in  f_anno.read().strip().split('\n')]
    
# 筛选所需要的的基因位置
tx_list = []
for line_anno_list in anno_list:
    if line_anno_list[-1] in target_gene_list:
        tx_list.append(line_anno_list[:3] + line_anno_list[-1:])
# 写出筛选后的 bed&gene 对应文件
with open(out_file, 'w') as f_out:
    for out_line_list in tx_list:
        f_out.write('\t'.join(out_line_list) + '\n')

调整后的格式样式：

2.bedtools 进行注释

代码如下（示例）：

bedtools intersect -a 44.bed -b gs.anno.tx -loj > out_file_loj
awk -v OFS="\t" "{print $1,$2,$3,$8}" out_file_loj | uniq > out_file_loj_uniq
bedtools groupby -i out_file_loj_uniq -g 1,2,3 -c 4 -o collapse > out_file_loj_uniq_collapse

结果样式：

总结

RefGene下载地址：http://genome.ucsc.edu/cgi-bin/hgTables