利用bedtools预测chip_seq数据的靶基因

欢迎关注”生信修炼手册”!

通常在分析peak区域对应的靶基因时，会选取转录起始位点TSS上下游一定长度的区域作为候选的靶基因范围，本文介绍下如何利用bedtools来对peak与TSS区域的overlap情况进行分析，从而得到靶基因，可以分为以下几步

1. 得到物种对应的TSS位点信息

以​​hg38​​​为例，通过UCSC的FTP服务可以得到物种对应的​​refFlat​​文件，链接如下

​​http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/​​

​​refFLat​​​和​​refGene​​​这两个文件记录的信息相同，​​refFlat​​​文件列数更少，这里我们选择下载​​refFlat.txt.gz​​, 该文件的内容如下所示

在原始文件中是没有第一行的标题的，我手动添加的标题是为了方便描述每列的含义，从该文件中可以得到TSS位点信息。

2. 整理TSS位点信息

bedtools要求输入的文件格式为​​bed​​​, ​​gff​​​, ​​vcf​​​等，这里我们需要把上述下载的原始文件转换为​​bed​​格式，用法如下

awk '{print $3"\t"$5"\t"$5"\t"$2"\t"$1"\t"$4}' > hg38.tss.bed

内容如下所示

3. 运行bedtools window

bedtools windows和intersect的功能类似，都是用于求两个区间A和B的交集，只不过​​window​​会在A区间的上下游加上一个可以自定义的长度之后，再与B区间求交集，原理示意如下

以TSS上下游5kb为例，用法如下

bedtools window -a hg39.tss.bed -b peak.bed -w 5000 -sm > overlap.txt

通过​​window​​这个命令，可以灵活的定义TSS上下游的区间，快速得到peak对应的靶基因。

·end·

—如果喜欢，快分享给你的朋友们吧—

扫描关注微信号，更多精彩内容等着你！