一、概述

1、概述

BEDTools是可用于genomic features的比较，相关操作及进行注释的工具。而genomic features通常使用Browser Extensible Data (BED) 或者 General Feature Format (GFF)文件表示，用UCSC Genome Browser进行可视化比较。bedtools总共有二三十个工具/命令来处理基因组数据。比较典型而且常用的功能举例如下：格式转换，bam转bed（bamToBed），bed转其他格式（bedToBam，bedToIgv）；对基因组坐标的逻辑运算，包括：交集（intersectBed，windowBed），”邻集“（closestBed），补集（complementBed），并集（mergeBed），差集（subtractBed）;计算覆盖度（coverage）（coverageBed，genomeCoverageBed）；
此外，还有一些强大而实用的工具（shuffleBed，groupBy，annotateBed，……）

2、与BEDTools使用相关的基本概念

已有的一些genome features信息一般由BED格式或者GFF格式进行存储。

1）genome features: 功能元素（gene），遗传多态(SNPs,INDELs or structural variants), 已经由测序或者其他方法得到的注释信息，也可以是自定义的一些特征信息。

2）genome features的基本信息：染色体或者scaffold的位置，起始位置，终止位置，哪条链，feature的name

3）Overlapping/intersecting features:两个genome features的区域至少有一个bp的共同片段

3、BED和GFF文件的一个差异

BED文件中起始坐标为0，结束坐标至少是1；GFF中起始坐标是1而结束坐标至少是1。

4、相关格式

1) BED format

BEDTools主要使用BED格式的前三列，BED可以最多有12列。BED格式的常用列描述如下：

chrom: 染色体信息， 如chr1, III, myCHrom, contig1112.23, 必须有

start: genome feature的起始位点，从0开始， 必须有

end: genome feature的终止位点，至少为1， 必须有

name: genome feature的官方名称或者自定义的一个名字

score: 可以是p值等等一些可以刻量化的数值信息

strands: 正反链信息

2) BEDPE format

可以用于描述不连续的genome features， 例如structural variations或者paired-end sequence alignments。和BED文件格式相比，就是一个记录中要有两个chrom, start，end

3) GFF format

类似于BED

4) genome files

BEDTools中的一些工具（genomeCoverageBed, complementBed, slopBed）需要物种的染色体大小的信息，genome file一般就是每行都是tab隔开，两列，一列为染色体的名字，第二列为这个染色体的大小。一般常用物种的genome file在BEDTools安装目录的/genome里面。

5) SAM/BAM format

BEDTools的两个工具：intersectBed，pairToBed支持BAM格式的输入和输出。有两个工具有助于：

Find BAM alignments that overlap (or not) with BED annotation and report them in BED format

Create a new BAM file of BAM alignments that overlap (or not) with BED annotations. This serves as a powerful way to refine alignment datasets based on biological interest.

二、BEDTools suite使用详细

1、intersectBed

用来求两个BED或者BAM文件中的overlap，overlap可以进行自定义是整个genome features的overlap还是局部。 加-wa参数可以报告出原始的在A文件中的feature，加-wb参数可以报告出原始的在B文件中的feature, 加-c参数可以报告出两个文件中的overlap的feature的数量，参数-s可以得到忽略strand的overlap，具体案例如下：

案例一：包含着染色体位置的两个文件，分别记为A文件和B文件。分别来自于不同文件的染色体位置的交集是什么？

$ cat A.bed
chr1 10 20
chr1 30 40
$ cat B.bed
chr1 15 25
$ bedtools intersect -a A.bed -b B.bed 
chr1 15 20 

案例二：包含着染色体位置的两个文件，分别记为A文件和B文件。求A文件中哪些染色体位置是与文件B中的染色体位置有overlap.
$ cat A.bed
 chr1 10 20
 chr1 30 40
 $ cat B.bed
 chr1 15 25
 $ bedtools intersect -a A.bed -b B.bed -wa
 chr1 10 20 

案例三：包含着染色体位置的两个文件，分别记为A文件和B文件。求A文件中染色体位置与文件B中染色体位置的交集，以及对应的文件B中的染色体位置.
 $ cat A.bed
 chr1 10 20
 chr1 30 40
 $ cat B.bed
 chr1 15 25
 $ bedtools intersect -a A.bed -b B.bed -wb
  chr1 15 20 chr1 15 25

案例四（经用）： 包含着染色体位置的两个文件，分别记为A文件和B文件。求对于A文件的染色体位置是否与文件B中的染色体位置有交集。如果有交集，分别输入A文件的染色体位置和B文件的染色体位置；如果没有交集，输入A文件的染色体位置并以'. -1 -1'补齐文件。

$ cat A.bed

chr1 10 20
chr1 30 40
$ cat B.bed
chr1 15 25
$ bedtools intersect -a A.bed -b B.bed -loj
chr1 10 20 chr1 15 25
chr1 30 40 . -1 -1

案例五：  包含着染色体位置的两个文件，分别记为A文件和B文件。对于A文件中染色体位置，如果和B文件中染色体位置有overlap,则输出在A文件中染色体位置和在B文件中染色体位置，以及overlap的长度.
 $ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed
chr1 15 20
chr1 18 25
$ bedtools intersect -a A.bed -b B.bed -wo
chr1 10 20 chr1 15 20 5
chr1 10 20 chr1 18 25 2 

案例六：  包含着染色体位置的两个文件，分别记为A文件和B文件。对于A文件中染色体位置，如果和B文件中染色体位置有overlap,则输出在A文件中染色体位置和在B文件中染色体位置，以及overlap的长度；如果和B文件中染色体位置都没有overlap,则用'. -1-1'补齐文件

$ cat A.bed
chr1 10 20
chr1 30 40
$ cat B.bed
chr1 15 20
chr1 18 25
$ bedtools intersect -a A.bed -b B.bed -wao
chr1 10 20 chr1 15 20 5
chr1 10 20 chr1 18 25 2 
chr1 30 40 . -1 -1

 案例七：  包含着染色体位置的两个文件，分别记为A文件和B文件。对于A文件中染色体位置，输出在A文件中染色体位置和有多少B文件染色体位置与之有overlap.

$ cat A.bed
 chr1 10 20
 chr1 30 40
 $ cat B.bed
 chr1 15 20
 chr1 18 25
 $ bedtools intersect -a A.bed -b B.bed -c
 chr1 10 20 2
 chr1 30 40 0 

 案例八(常用)：  包含着染色体位置的两个文件，分别记为A文件和B文件。对于A文件中染色体位置，输出在A文件中染色体位置和与B文件染色体位置至少有X%的overlap的记录。

$ cat A.bed
chr1 100 200
$ cat B.bed
chr1 130 201
chr1 180 220
$ bedtools intersect -a A.bed -b B.bed -f 0.50 -wa -wb  
chr1 100 200 chr1 130 201 

2、pairToBed

将BED PE或者paired-end的BAM文件与BED文件进行比较，搜索overlaps。默认结果如下图

       加-type both报告A中两个end都在B中存在overlap的；

加-type notboth报告A中任一end都在B中存在overlap的； 

加-type ispan报告A的pair是否是横跨B的两侧，

还有类似的-type ospan, -type notispan；

加参数-f可以指定最小的overlap的比例，超过这个比例则报告；

参数-s可以得到忽略strand的

举例：

overlappairToBed -a A.bedpe -b B.bed -type bothpairToBed -a A.bedpe

 -b B.bed -f 0.5pairToBed -abam pairedReads.bam -b simreps.bed -bedpe

3、bedtools multicov.

利用比对的bam文件对每个基因进行read数目计数。

首先注意它需要sort的bam文件及bam的index

bedtools multicov [OPTIONS] -bams BAM1 BAM2 BAM3 … BAMn -bed  <BED/GFF/VCF>

例子：

bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed interest.bed

interest.bed这个文件是必须的，可能需要自己制作，其实用gtf文件也可以的

chr1 0   10000   gene1

chr1 10000   20000   gene2

chr1 20000   30000   gene3

chr1 30000   40000   gene4

输出结果前三列是坐标，第四列是基因名，跟我们的bed文件一样，只是最后三列是三个样本的计数，是添加上来的！

chr1 0       10000   gene1    100 2234    0

chr1 10000   20000  gene2    123 3245    1000

chr1 20000   30000   gene3    213 2332    2034

chr1 30000   40000   gene4    335 7654    0

 得到的这个矩阵就可以去用DESeq包来进行差异分析啦！

4.bedtools merge

 用于合并位于同一个bed/gff/vcf 文件中的重叠区域。

 Bedtools merge  [OPTION] –i <bed/gff/vcf>

-s 必须相同(正负)链的区域才合并（软件默认不考虑正负链特征）

-n 报告合并的区域数量，没有被合并则1

-d 两个独立区域间距小于（等于）该值时将被合并为一个区域

-nms 报告被合并区域的名称

-scores 报告几个被合并特征区域的scores

5、其他

1）pairToPair

比较BEDPE文件搜索overlaps, 类似于pairToBed。

2）bamToBed

将BAM文件转换为BED文件或者BEDPE文件。

bamToBed -i reads.bam

3）windowBed类似于intersectBed, 但是可以指定一个数字，让A中的genome feature增加上下游去和B中的genome features进行overlap。默认情况这个值为1000，可以使用-w加定义，可以用-l指定是上游，用-r指定下游

windowBed -a A.bed -b B.bed -w 5000

windowBed -a A.bed -b B.bed -l 200 -r 20000

4）subtractBed在A中去除掉B中有的genome features

5）coverageBed可以计算深度和覆盖度。如计算基因组任意1Kb的测序read的覆盖度

6）genomeCoverageBed。可以计算给定bam文件在基因组上的覆盖度及每个碱基的深度。

三、软件相关论文：

Quinlan, A.R. & Hall, I.M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841-842 (2010).