知识中心 - 北京概普生物科技有限公司(GapTech)

一文掌握bed运算

生信干货青壬 ·2018年3月14日 02:01

1. bed格式概述

bed是二代测序数据后期常见的数据格式，用来表述基因组区段位置及其附属信息。最基本的BED3格式即：chrom - start - end。除前三列外的列信息均为附属信息，形成拓展的bed格式。可增加的附属信息最多有9列，在UCSC官网上有相关定义如下：一般来说常见的bed格式有以下几种：

bed	format	example
bed3	chrom, start, end	chr4 11573 14409
bed4	chrom, start, end, name	chr4 11573 14409 uc001aaa.3
bed5	chrom, start, end, name, score	chr4 11573 14409 uc001aaa.3 0
bed6	chrom, start, end, name, score,strand	chr4 11573 14409 uc001aaa.3 0 +
bed12	chrom, start, end, name, score,strand	chr4 11573 14409 uc001aaa.3 0 + 11873 11873 0 3 354,109,1189, 0,739,1347,

除此之外，还有很多基于bed的变异格式，如bedGraph等。

关于bed的一些冷知识：

bed是0起序的(zero-based)，即chrom1的第一个碱基写作 chr1 0 1；
bed的染色体名称并不作统一命名要求，需要分析者自己规范；
bed默认采用tab键分割(tab-delimited)。

2. bedtools的运用

在生信发展的早期，对于bed格式文件的处理需要分析者自己进行编程以进行坐标运算。其思路比较简单，即利用线性的扫描不断移动指针最终对一个个染色体完成运算。目前主流使用的bed运算工具是由犹他大学Quinlan组开发的bedtools，该工具操作简便且运行效率高，逐渐成为二代测序数据分析者的瑞士军刀。以下将对bedtools的用法进行梳理介绍：

head cpg.bed 
chr1    28735   29810   CpG:_116
chr1    135124  135563  CpG:_30
chr1    327790  328229  CpG:_29
chr1    437151  438164  CpG:_84
chr1    449273  450544  CpG:_99
head exons.bed 
chr1    11873   12227   NR_046018_exon_0_0_chr1_11874_f 0   +
chr1    12612   12721   NR_046018_exon_1_0_chr1_12613_f 0   +
chr1    13220   14409   NR_046018_exon_2_0_chr1_13221_f 0   +
chr1    14361   14829   NR_024540_exon_0_0_chr1_14362_r 0   -
chr1    14969   15038   NR_024540_exon_1_0_chr1_14970_r 0   -

1. bed做交集运算

bed的运算本质上是集合运算，因此基于数学性质继承下来的交并补是必不可少的。

相信凡是使用过bedtools的同学都对这张调皮的手绘风格示意图印象深刻。这张图直白展示了bedtools交集运算的参数作用。

最基础的交集运算如下所示：

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

-wa -wb 可以对每一次做交集运算附带上相应bed文件。注：一般来说这个参数的设置是为了获取附加列信息与计算出的交集的关系。

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        -
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   +

-v 相当于集合运算的A-B，即不与B重叠的A部分。 -wo 会对重叠长度进行统计

bedtools intersect -a cpg.bed -b exons.bed -wo \
| head -2
chr1    28735   29810   CpG:_116    chr1    29320   29370   NR_024540_exon_10_0_chr1_29321_r        -   50
chr1    135124  135563  CpG:_30 chr1    134772  139696  NR_039983_exon_0_0_chr1_134773_r    0       439

-c 则可以对A包含多少个B进行统计

bedtools intersect -a cpg.bed -b exons.bed -c \
| head -3
chr1    28735   29810   CpG:_116    1
chr1    135124  135563  CpG:_30 1
chr1    327790  328229  CpG:_29 3

一般来说所谓重叠的标准低至1bp，我们可以通过-f参数进行调整。注意：-f表示的是重叠占A的比例，而不是绝对的长度。

bedtools intersect -a cpg.bed -b exons.bed \
-wo -f 0.50 \
| head -2
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

2. bed做并集运算

对于并集运算，比较tricky的一点在于拼接距离。如图所示，通过-d参数可以设定此阈值。

bedtools merge -i exons.bed -d 1000 -c 1 -o count | head -2
chr1    11873   18366   
chr1    24737   24891

我们还可以对指定的bed在交集中的个数进行统计：

bedtools merge -i exons.bed -c 1 -o count | head -n 2
chr1    11873   12227   1
chr1    12612   12721   1

除此之外，-o还提供了一系列聚合函数，我们可以通过manual详细参看。例如collapse：

bedtools merge -i exons.bed -d 90 -c 1,4 -o count,collapse | head -3
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

3. bed做补集运算

补集运算相对比较简单。相对于intersect -v来说，补集的计算背景是全集。因此我们需要提供计算基因组大小的txt文件，例如：

head genome.txt
chr1    249250621
chr10   135534747
chr11   135006516
chr12   133851895
chr13   115169878
chr14   107349540
chr15   102531392
chr16   90354753

基于此进行补集运算：

bedtools complement -i exons.bed -g genome.txt | head -n 4 
chr1    0   11873
chr1    12227   12612
chr1    12721   13220
chr1    14829   14969

4. 运用bedtools计算不同区域的genome coverage

直接计算出各个染色体的genome coverage频率分布

bedtools genomecov -i exons.bed -g genome.txt | head -5
chr1    0   241996316   249250621   0.970896
chr1    1   4276763 249250621   0.0171585
chr1    2   1475526 249250621   0.00591985
chr1    3   710135  249250621   0.00284908
chr1    4   388193  249250621   0.00155744

最终应当注意的是：bedtools的输入必须经过sort，并且需要统一各个bed染色体名称的表述（chr，chrom等）。通常利用：

sort -K1V -K2n

就能很好的完成排序操作。

欢迎关注生信人

点击以下「关键词」，查看往期内容：

一个物种一个家

TCGA | 小工具 | 数据库 |组装| 注释 | 基因家族 | Pvalue

基因预测 |bestorf | sci | NAR | 在线工具 | 生存分析 | 热图

舞台|基因组 | 黄金测序 | 套路 | 杂谈组装 | 进化 | 测序简史