用bedtools和bedops对bed文件进行各种运算
我经常和别人说 50% 高通量数据的后期分析概括一下都是各种位置之间的统计计算，因为但凡分析有参数据就会把所有信息都铆钉到参考基因组这个坐标系中，如果没有基因组，那就拼一个？这也是为什么 bedtools 从一推出就如此受欢迎。在生物信息领域，一个软件是不是会被大量使用，软件是不是真写的又快又好又规范都不是最重要的，最重要的是它是不是真的有用，然后再加上一点点出现的时机和运气。
在平时工作中，很多人会碰到类似于下面的场景： 已有两组感兴趣的基因组对应位置，这些位置可能是某些基因所代表的坐标，可能是某些分析找出的各种结合位点，例如 ChIP-seq 或者 ATAC-seq，还可能是一些基因组上的变异信息比如 indel 或者 snp。有了这两组位置的list，如果想知道一些交并补集的逻辑运算比较简单，在 bedtool 中常规的 intersect, closest, merge 等就可以完成。
如果想对一组位置在另一组中的某些数据进行相对复杂些的计算又要如何操作呢？比如你想知道在1024个上调基因位置区间内每个基因里有多少个某转录因子的peak位置，这些 peak 的最大值是多少，亦或者这些 peak 峰度的均值和变异系数是什么。
可以使用 bedtool groupby
在 bedtools 中有一个对应的命令叫做 groupby,官网文档对这个命令的基本解释如下：
Given a file or stream that is sorted by the appropriate “grouping columns” (-g), groupby will compute summary statistics on another column (-c) in the file or stream. This will work with output from all BEDTools as well as any other tab-delimited file or stream.
groupby 主要参数为 -g、-c 和 -o。其中，-g 指定根据那些信息进行group，默认是根据 bed 文件的 1-3 列； -c 是指定针对哪一列以 group 为单位进行各种运算； -o 则是指定进行何种运算。
目前，groupby 可以记性的运算内容如下：
sum - numeric only
count - numeric or text
count_distinct - numeric or text
min - numeric only
max - numeric only
mean - numeric only
median - numeric only
mode - numeric or text
antimode - numeric or text
stdev - numeric only
sstdev - numeric only
collapse (i.e., print a comma separated list) - numeric or text
distinct (i.e., print a comma separated list) - numeric or text
concat (i.e., print a comma separated list) - numeric or text
freqasc - print a comma separated list of values observed and the number of times they were observed.
Reported in ascending order of frequency*
freqdesc - print a comma separated list of values observed and the number of times they were observed.
Reported in descending order of frequency*
first - numeric or text
last - numeric or text
具体的使用示例去参考官方文档吧，懒得在这里浪费时间。只提一句，一般 groupby 和 intersect 连用比较常见。
推荐使用 bedops bedmap
在 bed 格式文件处理领域，bedtools 可谓一家独大。其实还有一个很强大的工具并没有引起太多人的注意，这个软件叫做bedops。bedops 最突出的两个特点，一个是快且省内存，一个是很灵活。他的速度要比之前版本的bedtools 快非常多（据说后来 bedtools 升级速度有所改善）同时很省内存，而且支持各种 pipe 的结合连用。另外一个隐藏优点是 bedtools 在处理不少大基因组文件时会因为超过软件本身的各种限制而报错无法使用，但 bedops 则不存在这类问题。
bedops 这款优秀的软件这里不做太多介绍，具体可以了解官方文档。现在只要知道 bedtools 中 95% 的工具在 bedops 中都有相应很方便的方法，而且 bedops 还有一些其他好用的命令。
针对本文提到的需求，在 bedops 中有一个对应的核心工具叫做bedmap，它可以根据你的基因组输入文件提供各种各样的统计操作。
bedmap 的输入文件为两个bed文件，其中一个用作分 group 的 ref-bed 文件可以是最简单的三列格式用来提供位置信息，另一个需要进行后续计算的文件 map-bed 则根据实际计算需求至少为4列或者5列，其中 1-3 列为位置信息，第 4 列为 id，第 5 列为要计算和统计的值。
bedmap 参数整体可以分为三大类。一类是处理参数，例如指定处理速度和数据格式等；第二类是 overlap 参数，即满足怎样的条件才进行计算，例如指定位置间 overlap 的碱基数或者比例等；第三类参数则为运算参数，其中一个子类是分值类计算，例如最大值，均值，中位数和变异系数等等，一个子类是非分值类参数，例如包含个数等。需要提醒的是，如果进行分值类运算，map bed 需要有5列，如果是非分值类运算 map bed 文件根据要计算的内容不同还略有差别，详情和使用示例可以参考官方的文档说明。
写在最后
就个人使用体验来说，我目前已经把绝大部分 bed 文件相关处理工作从 bedtools 转移到了 bedops，现在也把这个工具推荐给你。
当然，除了上面两个工具外，在 Python 和 R 中也有人进行了相应移植工作，以满足不同的使用环境和需求，在文章的最后列出来供你参考。
- Python：pybedtools; bedtools-python
- R: bedr
本文作者：思考问题的熊
版权声明：本博客所有文章除特别声明外，均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。
