2019 和豆豆一起跟着官网学习 bedtools

豆豆写于19.12.6
学了一天，终于学完了官网的教程，细细品味

我承认，我饿了

最近在学习ChIPseq分析，这个工具对数据格式转换和操作基因组坐标是很有帮助的。不过虽然很受欢迎，但还是有很多类似的工具可以用，比如bedops。慢慢来，不着急，先学下这个目前引用超过8500的”瑞士军刀“

0-1 前言

这个工具的官网在：https://bedtools.readthedocs.io/en/latest/

作者 Aaron R. Quinlan 2010年把这个工具发在了Bioinformatics上，他当时还在美国弗吉尼亚大学医学院，现在工作在犹他大学【他的实验室链接：http://quinlanlab.org/，这个网站做的真心不错】

点开官网链接，它的logo早已说明一切。它支持对BAM、BED、GFF/GTF、VCF等和基因组区域有关文件的intersect, merge, count, complement, and shuffle等操作。或许这些名词还看不懂，没关系，一点点跟着例子走：

0 下面👇正文开始

官方教程在：http://quinlanlab.org/tutorials/bedtools/bedtools.html
我也是跟着这个学一遍，记录自己的心得

软件安装不多说，conda直接搞定

0 测试数据的下载

# 创建一个属于你自己的学习目录
mkdir -p ~/biosoft/bedtools-learn
cd ~/biosoft/bedtools-learn
# 然后下载测试数据，一个压缩包，一个txt，4个bed文件
curl -O https://s3.amazonaws.com/bedtools-tutorials/web/maurano.dnaseI.tgz
curl -O https://s3.amazonaws.com/bedtools-tutorials/web/cpg.bed
curl -O https://s3.amazonaws.com/bedtools-tutorials/web/exons.bed
curl -O https://s3.amazonaws.com/bedtools-tutorials/web/gwas.bed
curl -O https://s3.amazonaws.com/bedtools-tutorials/web/genome.txt
curl -O https://s3.amazonaws.com/bedtools-tutorials/web/hesc.chromHmm.bed
# 【当然，如果自己的网络不给力，也可以在公众号回复”bedtools“即可】

tar -zxvf maurano.dnaseI.tgz
rm maurano.dnaseI.tgz

0.1 这些数据都是什么？

• 可以看到解压后有20个以f开头的bed文件（f 代表 fetal tissue 胎儿组织），表示来自大脑、心脏、肠、肾、肺、肌肉、皮肤和胃的20个不同胎儿组织样本。

• cpg.bed表示人类基因组的CpG岛

• exons.bed表示来自RefSeq的人类基因外显子.

• gwas.bed表示用GWAS鉴定的与人类疾病相关的SNPs

• hesc.chromHmm.bed表示基于ENCODE的ChIP-seq实验，利用chromHMM预测的人类胚胎干细胞基因组区间的功能

这四个文件都是从：UCSC Table Browser 获得的，把它们导入IGV是这样：

0.2 bedtools概况

这个工具有几十个子命令，比如使用时要输入：

bedtools intersect
bedtools merge
bedtools subtract
...

查看版本：

bedtools --version

强烈建议下面的命令自己手动敲一遍！

1 第一个最常用的：bedtools intersect

它的作用是比较两个或多个BED/BAM/VCF/GFF，然后找它们重叠的区域（就是指至少存在1个公共的碱基）

1.1 官方给出的手绘图感觉最好理解

如果只比较两个文件 A和B

解释一下”绿色部分的结果“

• 第一种没有参数的：最好理解，就是简单的返回A和B共有的区域

• 第二种加了-wa参数：它的含义是Write the original entry in A for each overlap，也就是把原始A文件的坐标输出来

• 第三种加了-v参数：只输出A中有且B中没有的区域（很像grep -v）

如果比较多个文件 A 和 B1、B2

• 第一种结果：输出A与B1、B2三者共有的区域

• 第二种加了-wa、-wb参数：先输出A的坐标，然后输出与A有重叠的B1、B2的区域坐标

1.2 默认的输出结果

通过-a、-b指定文件A、B

bedtools intersect -a cpg.bed -b exons.bed | head -5

chr1    29320   29370   CpG:_116
chr1    135124  135563  CpG:_30
chr1    327790  328229  CpG:_29
chr1    327790  328229  CpG:_29
chr1    327790  328229  CpG:_29

需要注意的是：这里打印出来的结果并不是原始的CpG记录的区间，而是它和exon.bed的重叠区间，因此坐标只是原来CpG.bed的子集

1.3 想看具体重叠双方的区域

Not only showing you where the intersections occurred, it shows you what intersected

bedtools intersect -a cpg.bed -b exons.bed -wa -wb | head -5

chr1    28735   29810   CpG:_116    chr1    29320   29370   NR_024540_exon_10_0_chr1_29321_r    0   -
chr1    135124  135563  CpG:_30 chr1    134772  139696  NR_039983_exon_0_0_chr1_134773_r    0   -
chr1    327790  328229  CpG:_29 chr1    324438  328581  NR_028322_exon_2_0_chr1_324439_f    0   +
chr1    327790  328229  CpG:_29 chr1    324438  328581  NR_028325_exon_2_0_chr1_324439_f    0   +
chr1    327790  328229  CpG:_29 chr1    327035  328581  NR_028327_exon_3_0_chr1_327036_f    0   +

就第一行结果来说：

# 3.2中得到的是
chr1    29320   29370   CpG:_116
# 3.3 得到的是
chr1    28735   29810   CpG:_116    chr1    29320   29370   NR_024540_exon_10_0_chr1_29321_r    0   -

看到原始的CpG记录是28735 - 29810，只不过有一个外显子坐标（29320 - 29370）包含其中，于是3.2的结果只记录了29320 - 29370

1.4 想统计重叠的碱基数

使用参数-wo  （write overlap），会在每一行结果末尾统计这个重叠区域有多长

例如：下面👇结果结尾的50、439

bedtools intersect -a cpg.bed -b exons.bed -wo | head -5

chr1    28735   29810   CpG:_116    chr1    29320   29370   NR_024540_exon_10_0_chr1_29321_r    0   -   50
chr1    135124  135563  CpG:_30 chr1    134772  139696  NR_039983_exon_0_0_chr1_134773_r    0   -   439
chr1    327790  328229  CpG:_29 chr1    324438  328581  NR_028322_exon_2_0_chr1_324439_f    0   +   439
chr1    327790  328229  CpG:_29 chr1    324438  328581  NR_028325_exon_2_0_chr1_324439_f    0   +   439
chr1    327790  328229  CpG:_29 chr1    327035  328581  NR_028327_exon_3_0_chr1_327036_f    0   +   439

1.5 统计发生重叠的features数量

这个feature的含义很广，不限于基因，还可以是外显子、转录本等等与基因组有关的信息

使用参数-c ，就会统计对文件A中的每个feature，在文件B中有多少个features与之有交集

例如：下面👇结果结尾的1、1、3、0、0

bedtools intersect -a cpg.bed -b exons.bed -c | head -5

chr1    28735   29810   CpG:_116    1
chr1    135124  135563  CpG:_30 1
chr1    327790  328229  CpG:_29 3
chr1    437151  438164  CpG:_84 0
chr1    449273  450544  CpG:_99 0

1.6 统计没有重叠的区域

使用-v参数【这个应该很容易理解】

bedtools intersect -a cpg.bed -b exons.bed -v | head -5

chr1    437151  438164  CpG:_84
chr1    449273  450544  CpG:_99
chr1    533219  534114  CpG:_94
chr1    544738  546649  CpG:_171
chr1    801975  802338  CpG:_24

1.7 设定一个最小重叠比例

之前说过，软件默认只要二者有1bp交叉，也算重叠区域。但这个我们可以自行调整。

例如，设定至少满足文件A的百分之50要与文件B有重叠，才统计。使用参数-f （fraction）

bedtools intersect -a cpg.bed -b exons.bed -wo -f 0.5 | head -5

chr1    135124  135563  CpG:_30 chr1    134772  139696  NR_039983_exon_0_0_chr1_134773_r    0   -   439
chr1    327790  328229  CpG:_29 chr1    324438  328581  NR_028322_exon_2_0_chr1_324439_f    0   +   439
chr1    327790  328229  CpG:_29 chr1    324438  328581  NR_028325_exon_2_0_chr1_324439_f    0   +   439
chr1    327790  328229  CpG:_29 chr1    327035  328581  NR_028327_exon_3_0_chr1_327036_f    0   +   439
chr1    788863  789211  CpG:_28 chr1    788770  794826  NR_047519_exon_5_0_chr1_788771_f    0   +   348

然后以第一行结果为例，在IGV中查看，发现这个确实是满足的，因为文件A（cpg.bed）的100%都落在了文件B（exons.bed）中：

1.8 「小贴士」sort后的数据运行的更快哦😯

sort的含义是：对染色体和起始坐标进行排序（比如最简单的命令：sort -k1,1 -k2,2n；当然对bam文件排序可以用samtools的sort命令）

比如作者在这里做了一个演示：

time bedtools intersect -a gwas.bed -b hesc.chromHmm.bed > /dev/null
1.10s user 0.11s system 99% cpu 1.206 total

time bedtools intersect -a gwas.bed -b hesc.chromHmm.bed -sorted > /dev/null
0.36s user 0.01s system 99% cpu 0.368 total

而且文件越大，sort的优势越明显，也许这正是得到bam文件后sort的用意吧

1.9 一次寻找多个文件的重叠区域

在bedtools的2.21.0版本之后，便可以实现参数-b指定多个文件，大大简化了一次实验多个数据之间的比较

例如：寻找外显子与 CpG islands、GWAS SNPs、ChromHMM annotations的交叉【记得进行sort哦😯】：

bedtools intersect -a exons.bed -b cpg.bed gwas.bed hesc.chromHmm.bed -sorted | head -5

chr1    11873   11937   NR_046018_exon_0_0_chr1_11874_f 0   +
chr1    11937   12137   NR_046018_exon_0_0_chr1_11874_f 0   +
chr1    12137   12227   NR_046018_exon_0_0_chr1_11874_f 0   +
chr1    12612   12721   NR_046018_exon_1_0_chr1_12613_f 0   +
chr1    13220   14137   NR_046018_exon_2_0_chr1_13221_f 0   +

不过，这样我们看不到 文件（B1、B2…）与文件A到底在哪里产生了重叠

方案一：加上`-wa`、`-wb`参数

bedtools intersect -a exons.bed -b cpg.bed gwas.bed hesc.chromHmm.bed -sorted -wa -wb \
  | head -10000 \
  | tail -5

看到结果的第7列记录了一些数字：这代表了是文件B中的第几个文件（比如这个3 就表示hesc.chromHmm.bed这个文件，说明它和文件A exons.bed 存在交集）；另外，如果多个文件B与文件A的相同区域有交集，那么也是每一行记录一个文件

方案二：可能用数字表示文件还是不直观，我们还要回过头去查看到底是哪个文件

于是可以用参数-names提前给多个文件B起绰号

bedtools intersect -a exons.bed -b cpg.bed gwas.bed hesc.chromHmm.bed -sorted -wa -wb -names cpg gwas chromhmm \
  | head -10000 \
  | tail -10 

再看一下结果，这样是不是就方便查看了？

2 第二个常用的：bedtools merge

许多数据集的基因组feature坐标经常是连续的（比如ChIPseq的结果），就像下图的蓝色部分

于是可以把这些连续的基因组小区间连接起来，拼成一个连续的大区间

需要拼接的输入文件（bed/gff/vcf）必须是排序（sort）过的【不sort会报错🙅‍♂️】

先按染色体，再按起始位点，这样保证merge的算法执行起来非常顺畅，而基本不需要消耗内存再次加工

# 这样操作
sort -k1,1 -k2,2n foo.bed > foo.sort.bed

2.1 以exon.bed为例，展示merge的作用

head -n10 exons.bed

注意看：第3行和第4行的区间是有重叠的，因此它们可以进行merge

bedtools merge -i exons.bed | head -10

可以看到，merge之后原来的第3行（13220 - 14409）和第4行（14361 - 14829）坐标合并成了13220 - 14829

2.2 merge后统计之前重叠区间的数量

就像上面👆的第3行和第4行，它们合并在一起后，我们希望能直观看出新的区间13220 - 14829 是不是合并后的；如果是，合并了多少重叠的区域？

使用参数：-c和-o，其中-c指定按照第几列进行统计；-o指定计算方法（默认是sum）

例如：统计数量可以指定-o为count，此时按照第几列统计都是一样的【但如果指定-o sum，那么指定的列必须是数值，例如第2列 ”start position“，结果就会把重叠区域的start position 相加】

# 仅仅统计次数就用count【可以指定第一列chr，但实际上还是按第2-3列的坐标进行识别】
bedtools merge -i exons.bed -c 1 -o count | head -n 10

chr1    11873   12227   1
chr1    12612   12721   1
chr1    13220   14829   2
chr1    14969   15038   1
chr1    15795   15947   1
chr1    16606   16765   1
chr1    16857   17055   1
chr1    17232   17368   1
chr1    17605   17742   1
chr1    17914   18061   1

# 默认是sum【此时就要指定第2列或其他数值型的列】
bedtools merge -i exons.bed -c 2 -o sum | head -n 10

chr1    11873   12227   11873
chr1    12612   12721   12612
chr1    13220   14829   27581
chr1    14969   15038   14969
chr1    15795   15947   15795
chr1    16606   16765   16606
chr1    16857   17055   16857
chr1    17232   17368   17232
chr1    17605   17742   17605
chr1    17914   18061   17914

2.3 即使不重叠，也能merge

使用-d参数（distance），就像第4部分刚开始一样，指定了 -d 10，一开始相距10bp的两个区间也合并在了一起。

再举个例子：如果认为相聚不超过1000bp的两个区间，依然表示同一个，那么也可以合并

bedtools merge -i exons.bed -d 1000 -c 1 -o count | head -10

chr1    11873   18366   12
chr1    24737   24891   1
chr1    29320   29370   1
chr1    34610   36081   6
chr1    69090   70008   1
chr1    134772  140566  3
chr1    323891  328581  10
chr1    367658  368597  3
chr1    621095  622034  3
chr1    661138  665731  3

2.4 merge了这么多，我怎么知道是哪些？

还是借助上面提到的-o参数，它除了能统计个数、计算值的大小，还能输出名称【使用collapse】。只要给它指定了对哪一列进行操作

例如：exon.bed中的第4列是exon的名称，这样在count结束后，就可以再根据第4列输出名称【如果存在合并的现象，那么其中所有的exon都会列出来，并以逗号分隔】

bedtools merge -i exons.bed -d 90 -c 1,4 -o count,collapse | head -10

chr1    11873   12227   1   NR_046018_exon_0_0_chr1_11874_f
chr1    12612   12721   1   NR_046018_exon_1_0_chr1_12613_f
chr1    13220   14829   2   NR_046018_exon_2_0_chr1_13221_f,NR_024540_exon_0_0_chr1_14362_r
chr1    14969   15038   1   NR_024540_exon_1_0_chr1_14970_r
chr1    15795   15947   1   NR_024540_exon_2_0_chr1_15796_r
chr1    16606   16765   1   NR_024540_exon_3_0_chr1_16607_r
chr1    16857   17055   1   NR_024540_exon_4_0_chr1_16858_r
chr1    17232   17368   1   NR_024540_exon_5_0_chr1_17233_r
chr1    17605   17742   1   NR_024540_exon_6_0_chr1_17606_r
chr1    17914   18061   1   NR_024540_exon_7_0_chr1_17915_r

3 学习第三个：bedtools complement 实现反选

给定一个feature坐标信息文件，我们如果不关心其中标记的区间，而是想看看有哪些区间不在这个文件中

例如，有一个ChIP-seq的peaks信息，现在想知道有哪些区域没有被抗体结合，就可以用complement

现在有了外显子的bed文件，通过反选，我们就能获得内含子或基因间区的坐标。那么既然是反选，除了exon.bed，还要有一个总体的范围（也就是基因组各个染色体的长度信息）【用-g指定】

# 首先要确定整体的操作范围
head -10 genome.txt

chr1    249250621
chr10    135534747
chr11    135006516
chr11_gl000202_random    40103
chr12    133851895
chr13    115169878
chr14    107349540
chr15    102531392
chr16    90354753
chr17    81195210

# 实现反选，选出不属于exon的区域
bedtools complement -i exons.bed -g genome.txt | head -10

chr1    0   11873
chr1    12227   12612
chr1    12721   13220
chr1    14829   14969
chr1    15038   15795
chr1    15947   16606
chr1    16765   16857
chr1    17055   17232
chr1    17368   17605
chr1    17742   17914

看到Chr1的第一个外显子就是从11874开始

但是不知道你注意没有，前面exon.bed中记录的chr1的第一个外显子坐标是

chr1 11873 12227 1 NR_046018_exon_0_0_chr1_11874_f

并非是11874开始，这是因为：

坐标系统 ：

• 1-based coordinate system ：序列的第一个碱基设为数字1，如SAM,VCF,GFF,wiggle格式

• 0-based coordinate system ：序列的第一个碱基设为数字0，如BAM, BCFv2, BED, PSL格式

6 学习第四个：bedtools genomecov 计算基因组覆盖度和深度

这个功能看上去摸不着头脑，实际上不难并且很实用

最简单的命令是：

bedtools genomecov -i exons.bed -g genome.txt

这个命令的默认输入是bed/gff/vcf ，并且需要加一个参考基因组的染色体大小记录文件。

6.1 来看一张关于genomecov 的图

看到这张图，第一反应是不是：这些0、1、2、3哪来的，有什么用？【没关系，一起探索】

6.2 什么是测序深度？

我们测序（尤其是ChIPseq），比较重要的一个目的就是查看测序深度，比如这些：

大家肯定对sam、bam文件不陌生，其中就记录了测序reads比对到基因组上的位置信息

【想对sam、bam文件补课？可以看：处理SAM、BAM你需要Samtools】

那么根据这个信息，就能反推基因组上每个碱基所对应的reads数，当然不用我们手动计算，samtools帮我们一把。可以看到下面得到的这个测序深度文件test.depth，第一列是染色体编号，第二列是染色体上碱基的位置，第三列是覆盖到这个位点的reads数目：

# 以果蝇的bam数据为例
samtools depth dm.sort.bam > test.depth

head -5 test.depth
2L    40  1
2L    41  1
2L    42  1
2L    43  1
2L    44  1

不过，虽然这个文件很易读，但每个碱基都记录于此，文件岂不是太大？（以人类为例，就需要记录30亿个碱基的信息，十分的冗余，因此得到的depth文件一般要比bam文件还要大几倍）

于是，为了体现测序深度的同时节省空间，又开发了两种新的文件格式：bedGraph和wiggle

6.3 bedGraph是什么？和depth文件有什么区别？

首先来看bedGraph文件的产生：

bedtools  genomecov -ibam dm.sort.bam -bg > test.bedgraph

如果要指定一个bam作为输入，就不需要使用-i和-g参数了。可以直接用-ibam指定输入bam文件；

然后-bg表示产生bedgraph文件

如果使用-bga 就表示：将feature之间的间隔区域记做窗口0（而-bg的结果中，这部分是不记录的）

看下它的内容

head test.bedgraph

2L    39  65  1
2L    65  66  2
2L    66  90  3
2L    90  116 2
2L    116 117 1
2L    255 256 1
2L    256 261 2
2L    261 284 3
2L    284 306 4
2L    306 307 3

区别

bedgraph文件将之前depth文件的每个碱基测序深度整合成为一个个的窗口，第四列记录的就是窗口中的测序深度，再看那个图，可以看到bedgraph将原来depth中具有相同测序深度的连续的碱基，合并成了一个小窗口：

最后，看一眼文件的大小差异

7 如果你认为bedtools只能做这些，那就不是瑞士军刀了

它的小功能应该有几十个，但最应该学习的是处理思想
官方给的建议是：Analytical power in bedtools comes from the ability to “chain” together multiple tools in order to construct rather sophisicated analyses with very little programming - you just need genome arithmetic!

下面将介绍不同的例子，来增加对它的认识

例1

给定一个`a.bam`和`b.regions.bed`，现在想看看region中的坐标在bam中没有体现的那部分

思路很清晰：先找bam中没有覆盖到的区域（用-bga 标记这部分为0），然后把0挑出来就是不在bam上的区域；然后intersect解决

bedtools genomecov -ibam a.bam -bga | awk '$4==0' | bedtools intersect -a regions -b - >result

例2

检查bam文件中比对到外显子的覆盖率

bedtools bamtobed -i reads.bam | \
bedtools coverage -a - -b exons.bed \
> exons.bed.coverage

再进一步， 查看properly-paired reads的覆盖度

这里的-uf 0x2 目的就是提取出来bam中的properly-paired reads，只不过用来二进制表示

关于flags可以看：https://davetang.org/muse/2014/03/06/understanding-bam-flags/

samtools view -uf 0x2 reads.bam | \
coverageBed -abam - -b exons.bed \
> exons.bed.proper.coverage

例3

对每条链分别计算覆盖度

# 先看+链（重点是根据转换的bed文件中的+/-那一列进行区分）
bedtools bamtobed -i reads.bam | \
grep \+$  | \
bedtools coverage -a - -b genes.bed \
> genes.bed.forward.coverage

# 然后看-链
bedtools bamtobed -i reads.bam | \
grep \-$ | \
bedtools coverage -a - -b genes.bed \
> genes.bed.reverse.coverage

例4

对自己记录的SNPs位点进行筛选，找到那些既不在dbSNP也不在1000 genomes中的

思路：核心就是对bed进行过滤（-v参数）

bedtools intersect -a snp.calls.bed -b dbSnp.bed -v | \
bedtools intersect -a - -b 1KG.bed -v | \
> snp.calls.novel.bed

例5

对外显子捕获效率的可视化

这个检测的标准可以参考：For example, covering 90% of the target region at 20X coverage may be one metric to assess your ability to reliably detect heterozygotes

思路：有一个测序reads比对的bam文件，有一个探针捕获区域的bed文件（比如：NimbleGen's V3 exome capture regions here ），然后用bedtools coverage命令，加-hist参数表示输出bed文件中各个feature在bam上覆盖度的直方图结果

直接看代码：https://www.gettinggeneticsdone.com/2014/03/visualize-coverage-exome-targeted-ngs-bedtools.html

例8

检测两个数据之间的相关性

引入一个新的bedtools工具jaccard，它会计算一个杰卡德相似性系数

结果是0.0 to 1. 0的值，数越小相关性越小

# 检测同一个样本的不同数据【系数是0.50637】
bedtools jaccard \
    -a fHeart-DS16621.hotspot.twopass.fdr0.05.merge.bed \
    -b fHeart-DS15839.hotspot.twopass.fdr0.05.merge.bed

intersection    union-intersection  jaccard n_intersections
81269248    160493950   0.50637 130852

# 再看不同的样本的不同数据【系数是0.170995】
bedtools jaccard \
    -a fHeart-DS16621.hotspot.twopass.fdr0.05.merge.bed \
    -b fSkin_fibro_bicep_R-DS19745.hg19.hotspot.twopass.fdr0.05.merge.bed

intersection    union-intersection  jaccard n_intersections
28076951    164197278   0.170995    73261

另外，还能分析更多的样本之间相关性

这个就看（官网http://quinlanlab.org/tutorials/bedtools/bedtools.html）翻到最底部

小结

叫小结一点都不夸张，上面的内容都是初级内容，是带领我们熟悉这个工具。就像搭积木，工具简单，但实现的方法千千万。

bedtools全部小功能列表

其中每一个小功能都有自己的链接，点进去，又是一片天地

https://bedtools.readthedocs.io/en/latest/content/bedtools-suite.html

bedtools的示范操作

帮你快速记忆，有点像cheatsheet

https://bedtools.readthedocs.io/en/latest/content/example-usage.html

常见的一些问题

很多都是安装相关的，既然用conda安装，这里就省去了很多麻烦

https://bedtools.readthedocs.io/en/latest/content/faq.html