最全Bedtools使用说明--只看本文就够了

2020-11-15 本文已影响0人生信小书童

1 Bedtools:一个强大的基因组算法工具集

Bedtools是由犹他大学昆兰实验室开发的基因组算法工具集，它堪称是基因组分析工具中的瑞士军刀。Bedtools可以对基因组广泛使用的数据格式BAM,BED,GFF/GTF,VCF进行处理，进行取交集、并集、补集、计数以及格式转变等操作。Bedtools被设计成许多独立的小工具，这些小工具自身只能处理简单任务，对于复杂的分析过程可以通过多个简单工具的组合进行处理。Bedtools也提供了详尽的教程，请参考
https://Bedtools.readthedocs.io/en/latest/index.html

2 bedtools备注

从2.28.0版开始，bedtools使用htslib库支持CRAM格式
除了BAM文件，bedtools默认所有的输入文件都以TAB键分割
除非使用-sorted选项，bedtools默认不支持大于512M的染色体
如果没有使用-sorted参数对染色体按编码顺序进行排序（e.g., sort -k1,1 -k2,2n ），则必须使用-g参数输入相同排序染色体
bedtools要求染色体命名方案在比较文件中是相同的（例如‘chr1’和‘1’不能同时存在）

3 Bedtools 须知（典范示例）：

基因组特征可以是功能元件（如基因），基因多态性（SNP、INDEL、SV）或者是由测序或者基因组浏览器发现或管理的其他注释，以及各个实验室或者研究者规定的传统注释。最基本的基因组特征是其所在的染色体，起始位置，终止位置，正负链特征。最广泛使用的基因组格式是BED文件和GFF文件。现有的许多物种的基因组注释可以很容易地从UCSC基因组浏览器的“表格浏览器”中以BED和GFF格式下载。

3.1 输入文件参数

下列例子中Feature A与Feature B有交集，但是与Feature C没有交集：

bedtools intersect –a snps.bed –b exons.bed

如果fileA是BAM格式，需要使用”-abam”选项

bedtools intersect –abam alignedReads.bam –b exons.bed

如果只需要输入一个基因组特征，则需要使用”-i”选项

bedtools merge –i repeats.bed

3.2 各种格式的起始位置

3.2.1 BED文件的起始位置以0开始，终止位置以1开始

具体来说，bedtools使用UCSC基因组浏览器的内部数据库约定，起始位置为0，结束位置为1，下面的BED特征代表的是1号染色体的一个碱基。这样存储特征的格式是为了方便计算特征的长度，如果起始位置以1开始，特征的长度((end-start)+1)，这样计算比计较复杂

chr1   0        1    first_base

3.2.2 GFF文件的起始和终止位置都是以1开始

Bedtools自身能设别GFF文件的正确位置，不需要再对GFF文件的初始位置减去1转变成bed格式

3.2.3 VCF文件的坐标是以1开始

和GFF文件类似，Bedtools自身能设别VCF文件的正确位置，不需要再对VCF文件的初始位置减去1转变成bed格式

3.3 大部分情况FileB文件会被加载到内存中

当比较两个文件的特征时，file B一般会被加入到内存中，file A会被遍历各行与内存中的file B进行比较。因此为了减少对内存的使用，一般需要将比较小的文件当做file B。例如，当比较一个含有数百万条reads的比对文件和数千条基因特征的注释文件时，会将比对文件当做file A，将注释文件当做file B。

3.4 Bedtools可以将特征以管道符的方式作为标准输入

当需要比较两个文件时，Bedtools将file A或者file B转化为”stdin” 或者 “-” 进行表示

cat snps.bed | bedtools intersect –a stdin –b exons.bed
cat snps.bed | bedtools intersect –a - –b exons.bed

当只需要输入一个文件时，Bedtools处理标准输入时可以忽略”-i”参数

cat snps.bed | bedtools sort –i stdin
cat snps.bed | bedtools sort

bedtools的结果直接以标准输出的方式输出

bedtools intersect –a snps.bed –b exons.bed
chr1 100100 100101 rs233454
chr1 200100 200101 rs446788
chr1 300100 300101 rs645678

以”>”符号输出到特定文件

bedtools intersect –a snps.bed –b exons.bed > snps.in.exons.bed
cat snps.in.exons.bed
chr1 100100 100101 rs233454
chr1 200100 200101 rs446788
chr1 300100 300101 rs645678

3.5 用”-h”参数来获取bedtools的帮助

使用“-h”参数，bedtools会列出使用的示例和全部的参数文件

3.6 BED文件位置信息不能包含负值，起始位置要<=终止位置

Bedtools通常不能包含负值。但是，在特殊情况下可以将BEDPE位置设置为-1,以指示BEDPE的一个或多个端点位置未对齐。

3.7 写出未压缩的BAM文件

当用一组复杂的管道处理大的BAM文件时，向下一个流程传递未压缩的BAM文件有利于节省解压的时间，所以用Bedtools输出BAM文件（intersect,window）时可以用-ubam选择性的输出未压缩的BAM文件

4 Bedtools各功能参数列表：

Utility	Description
annotate	注释多个文件的覆盖特征
bamtobed	将bam文件转换为BED格式
bamtofastq	将bam文件转换为FASTQ格式
bed12tobed6	将BED12间隔转换为BED6间隔
bedpetobam	将BEDPE间隔转化为BAM格式
bedtobam	将间隔转换为BAM记录
closest	寻找最近的潜在的非重叠区间
cluster	聚类（不合并）重叠的区间
complement	获得区间的补集
coverge	计算特定区间的覆盖
expand	根据列值重复行数
flank	从当前区间侧翼创建新的区间
genomecov	从整个基因组计算覆盖深度
getfasta	根据区间从FASTA文件中提取序列
groupby	按特征列进行分组
igv	创建一个IGV快照批处理脚本
intersect	用各种不同的方式寻找重叠区域
jaccard	计算b/w两个间隔区域的Jaccard统计量
links	创建一个连接UCSC位点的HTML页面
makewindows	创建基因组区间窗口
map	为每个重叠区间队列应用一个函数
maskfasta	利用区间隐藏FASTA文件序列
merge	合并重叠区间形成一个新的区间
multicov	计算多个BAM文件在特定区间的覆盖深度
multiinter	标记多个区间文件的公共区间
nuc	分析FASTA文件某区间的核酸含量
overlap	计算两个区间重叠范围的长度
pairtobed	找出以各种方式重叠区间的对
pairtopair	找出以各种方式重叠其他配对的配对
random	产生一个基因组的随机区间
reldist	计算两个文件的相对距离分布
shift	调整区间的位置
shuffle	在基因组中随机重新分配时间间隔
slop	调整区间的大小
sort	对区间进行排序
subtract	对两个区间文件取差集
tag	根据区间的重叠区域标记BAM Tag
unionbedg	根据多个BEDGRAPH文件合并覆盖区间
window	在一个间隔周围的窗口寻找重叠区间

4.1 annotate

<span id="mycode">Bedtools annotate</span> 命令用于统计一个BED/VCF/GFF文件与多个BED/VCF/GFF文件的重叠区域的比例和数目，这样就可以知道一个特征与多个特征的一致程度。
示例：

bedtools annotate [OPTIONS] -i <BED/GFF/VCF> -files FILE1 FILE2 FILE3 ... FILEn
$ cat variants.bed
chr1 100  200   nasty 1  -
chr2 500  1000  ugly  2  +
chr3 1000 5000  big   3  -

$ cat genes.bed
chr1 150  200   geneA 1  +
chr1 175  250   geneB 2  +
chr3 0    10000 geneC 3  -

$ cat conserve.bed
chr1 0    10000 cons1 1  +
chr2 700  10000 cons2 2  -
chr3 4000 10000 cons3 3  +

$ cat known_var.bed
chr1 0    120   known1   -
chr1 150  160   known2   -
chr2 0    10000 known3   +

$ bedtools annotate -i variants.bed -files genes.bed conserve.bed known_var.bed
chr1  100     200     nasty   1       -       0.500000        1.000000        0.300000
chr2  500     1000    ugly    2       +       0.000000        0.600000        1.000000
chr3  1000    5000    big     3       -       1.000000        0.250000        0.000000

参数：

Option	Description
-counts	输出的是多个特征与–i输入特征重叠的个数. 默认输出的是多个特征与-i输入特征重叠的比例
-both	输出的不仅是多个特征与–i输入特征重叠的个数，还包括多个特征与-i输入特征重叠的比例
-s	规定了相同的正负链。只有A和B首先具有相同的正负链时，才认为具有重叠区域，而默认不考虑正负链信息
-S	规定了相反的正负链. 只有A和B首先具有相反的正负链时，才认为具有重叠区域，而默认不考虑正负链信息

4.2 bamtobed

Bedtools bamtobed命令用于将序列比对文件bam格式转化为BED,BED12或者BEDPE格式。
示例：

bedtools bamtobed [OPTIONS] -i <BAM>
$ bedtools bamtobed -i reads.bam | head -3
chr7   118970079   118970129   TUPAC_0001:3:1:0:1452#0/1   37   -
chr7   118965072   118965122   TUPAC_0001:3:1:0:1452#0/2   37   +
chr11  46769934    46769984    TUPAC_0001:3:1:0:1472#0/1   37   -

参数：

Option	Description
-bedpe	将BAM格式转化为BEDPE格式。只有成对双端的reads才会被输出。当成对的双端reads分别比对到不同的染色体时，染色体字符低的染色体会被排在一行的前面。如果成对的reads,只有一条reads和参考序列比对上，那么没有比对上的那条reads的染色体和正负链会被设置为“.”，而起始和终止位置会被设置为-1。默认输出的格式是BED格式
-bed12	输出BED12格式。默认输出格式BED6格式
-tag	用BAM的其他tag当做BED score.默认的BED scroe值是bam文件中的MAPQ值
-cigar	输出cigar值当做BED文件的第7列

4.3 bamtofastq

bedtools bamtofastq 命令用于从比对的bam文件中提取fastq文件

示例：

bedtools bamtofastq [OPTIONS] -i <BAM> -fq <FASTQ>

参数：

Option	Description
-fq2	默认-fq输出的是一个fastq文件，添加-fq2参数可以将成对的fastq文件分别输出到两个文件中。但是输入的bam文件需要先对reads按名字进行排序 (samtools sort -n aln.bam aln.qsort)

4.4 bed12tobed6

bedtools bed12tobed6命令用于将BED12格式按照基因的外显子区块分隔成BED6格式

示例：

bedtools bed12tobed6 [OPTIONS] -i <BED12>
head data/knownGene.hg18.chr21.bed | tail -n 3
chr21 10079666  10120808   uc002yiv.1  0  -  10081686  1 0 1 2 0 6 0 8  0     4   528,91,101,215, 0,1930,39750,40927,
chr21 10080031  10081687   uc002yiw.1  0  -  10080031  1 0 0 8 0 0 3 1  0     2   200,91,    0,1565,
chr21 10081660  10120796   uc002yix.2  0  -  10081660  1 0 0 8 1 6 6 0  0     3   27,101,223,0,37756,38913,

head data/knownGene.hg18.chr21.bed | tail -n 3 | bed12ToBed6 -i stdin
chr21 10079666  10080194  uc002yiv.1 0  -
chr21 10081596  10081687  uc002yiv.1 0  -
chr21 10119416  10119517  uc002yiv.1 0  -
chr21 10120593  10120808  uc002yiv.1 0  -
chr21 10080031  10080231  uc002yiw.1 0  -
chr21 10081596  10081687  uc002yiw.1 0  -
chr21 10081660  10081687  uc002yix.2 0  -
chr21 10119416  10119517  uc002yix.2 0  -
chr21 10120573  10120796  uc002yix.2 0  -

参数：

Option	Description
-i	输入的是BED12格式. 通过管道输入时以”stdin”代替

4.5 bedtobam

bamtools bedtobam命令用于将bed格式转变成BAM格式。这对于储存大型压缩注释文件和可视化非常有用

示例：

bedtools bedtobam [OPTIONS] -i <BED/GFF/VCF> -g <GENOME> > <BAM>
head -5 rmsk.hg18.chr21.bed
chr21 9719768  9721892  ALR/Alpha  1004  +
chr21 9721905  9725582  ALR/Alpha  1010  +
chr21 9725582  9725977  L1PA3 3288 +
chr21 9726021  9729309  ALR/Alpha  1051  +
chr21 9729320  9729809  L1PA3 3897 -

bedToBam -i rmsk.hg18.chr21.bed -g human.hg18.genome > rmsk.hg18.chr21.bam

samtools view rmsk.hg18.chr21.bam | head -5
ALR/Alpha  0   chr21 9719769  255  2124M *  0  0  *  *
ALR/Alpha  0   chr21 9721906  255  3677M *  0  0  *  *
L1PA3      0   chr21 9725583  255  395M  *  0  0  *  *
ALR/Alpha  0   chr21 9726022  255  3288M *  0  0  *  *
L1PA3      16  chr21 9729321  255  489M  *  0  0  *  *

参数：

Option	Description
-mapq	设定输出BAM文件的mapping quality.默认的mapping quality是255
-ubam	输出未压缩的BAM文件，默认输出的是压缩后的BAM文件

4.6 closet

与intersect类似，closest命令会寻找A和B两个BED文件中的重合特征。但是，如果A和B的BED区域没有重叠，那么closet命令将会寻找B中与A的起始或者终止位置最近的区域。

示例：
bedtools closest [OPTIONS] -a <FILE>
-b <FILE1, FILE2, ..., FILEN>

cloest 示意图:

closest.png

4.7 Cluster

与merge功能类似，cluster命令会报告一个区间中的重叠区域。与merge命令不同，cluster命令不会合并重叠区间形成一个新区间，而是赋予各个区间一个聚类ID编号来反应各个区间的重叠或者远近。这对于反应一个区间中的重叠区域具有很好的效果。

示例：

bedtools cluster [OPTIONS] -i <BED/GFF/VCF>
$ cat A.bed
chr1  100  200
chr1  180  250
chr1  250  500
chr1  501  1000

$ bedtools cluster -i A.bed
chr1  100     200     1
chr1  180     250     1
chr1  250     500     1
chr1  501     1000    2

参数：

Option	Description
-s	具有相同的正负链才能归为一类.默认忽略正负链进行分类
-d	规定了两个特征聚类的最长距离. 默认的是有重叠区间或者两个区间紧邻时才能聚类

cluster 示意图：

cluster.png

4.8 Complement

bedtools complement 命令返回BED文件在一个基因组文件上的补集

示例：

bedtools complement -i <BED/GFF/VCF> -g <GENOME>
$ cat A.bed
chr1  100  200
chr1  400  500
chr1  500  800

$ cat my.genome
chr1  1000
chr2  800

$ bedtools complement -i A.bed -g my.genome
chr1  0   100
chr1  200 400
chr1  800 1000
chr2  0   800

complement 示意图：

complement.png

4.9 coverage

Bedtools coverage命令可以计算BED文件B在A文件上覆盖的深度和广度。例如bedtools coverage可以计算B文件在A文件特定的1kb窗口上的覆盖情况。Coverage命令的一大优势是它不仅计算与A文件各特征的重叠数目，而且计算重叠的比例。此外coverage也计算了与A文件各个区间重叠的区域的大小。

示例：

bedtools coverage [OPTIONS] -a <FILE> \
                             -b <FILE1, FILE2, ..., FILEN>
$ cat A.bed
chr1  0   100
chr1  100 200
chr2  0   100

$ cat B.bed
chr1  10  20
chr1  20  30
chr1  30  40
chr1  100 200

$ bedtools coverage -a A.bed -b B.bed
chr1  0   100  3  30  100 0.3000000
chr1  100 200  1  100 100 1.0000000
chr2  0   100  0  0   100 0.0000000

4.10 flank

bedtools flank 命令会在输入的区间的左右两侧形成两个新的区间，但是在两翼形成的区间不会超过染色体的大小范围（也就是说新形成的两翼区间起始位置要大于0，终止位置要小于染色体的大小）

示例：

bedtools flank [OPTIONS] -i <BED/GFF/VCF> -g <GENOME> [-b or (-l and -r)]
$ cat A.bed
chr1 100 200
chr1 500 600

$ cat my.genome
chr1 1000

$ bedtools flank -i A.bed -g my.genome -b 5
chr1  95      100
chr1  200     205
chr1  495     500
chr1  600     605

$ bedtools flank -i A.bed -g my.genome -l 2 -r 3
chr1  98      100
chr1  200     203
chr1  498     500
chr1  600     603

参数：

Option	Description
-b	两翼扩展范围的大小. 需是整数
-l	左翼扩展范围的大小 .需是整数
-r	右翼扩展范围的大小 .需是整数

4.11 genomecov

bedtools genomecov 用于在基因组层面计算各个位点的覆盖深度，覆盖深度以单个碱基展示（-d参数），覆盖深度以BED区域的形式展示（-bg参数）
需要注意的点是：

如果输入的格式是BED/GFF/VCF, 通过-i参数输入的文件必须对染色体和位置进行简单的排序（sort –k 1,1 in.bed），当输入的文件是BED/GFF/VCF格式时，必须提供-g genome文件
当输入的格式是BAM格式时，需要由-ibam参数导入，输入的BAM文件也需要先进行排序

示例：

bedtools genomecov [OPTIONS] [-i|-ibam] -g (iff. -i)
$ bedtools genomecov -ibam NA18152.bam -bg | head
chr1  554304  554309  5
chr1  554309  554313  6
chr1  554313  554314  1
chr1  554315  554316  6
chr1  554316  554317  5
chr1  554317  554318  1
chr1  554318  554319  2
chr1  554319  554321  6
chr1  554321  554323  1
chr1  554323  554334  7

genomecov 示意图：

genomecov.png

4.12 getfasta

Bedtools getfasta 命令可以从BED/GFF/VCF等格式定义的区间提取FASTA文件的中特定区域的序列。需要注意的点是：

输入的FASTA文件的表头必须和BED文件中染色体的列的名称精确匹配
可以用UNIX的fold –w 60命令规定每行显示的碱基的数目

示例：

bedtools getfasta [OPTIONS] -fi <input FASTA> -bed <BED/GFF/VCF>
$ cat test.fa
>chr1
AAAAAAAACCCCCCCCCCCCCGCTACTGGGGGGGGGGGGGGGGGG

$ cat test.bed
chr1 5 10

$ bedtools getfasta -fi test.fa -bed test.bed
>chr1:5-10
AAACC

# optionally write to an output file
$ bedtools getfasta -fi test.fa -bed test.bed -fo test.fa.out

$ cat test.fa.out
>chr1:5-10
AAACC

getfasta 示意图：

getfasta.png

4.13 groupby

Bedtools groupby 是一个有用的工具它可以对BED文件的部分列进行分组，并对指定的列进行统计。例如，如果一个文件的特定列（-g参数指定）已经进行了分组，groupby命令可以对指定的列（-c参数）进行统计。本工具不仅适用于基因组相关的研究，只要符合TAB键分隔的分组都可以用此工具处理。

示例：

bedtools groupby [OPTIONS] -i <input> -g <group columns> -c <op. column> -o <operation>
$ cat variantsToRepeats.bed
chr21  9719758 9729320 variant1   chr21  9719768 9721892 ALR/Alpha   1004  +
chr21  9719758 9729320 variant1   chr21  9721905 9725582 ALR/Alpha   1010  +
chr21  9719758 9729320 variant1   chr21  9725582 9725977 L1PA3       3288  +
chr21  9719758 9729320 variant1   chr21  9726021 9729309 ALR/Alpha   1051  +
chr21  9729310 9757478 variant2   chr21  9729320 9729809 L1PA3       3897  -
chr21  9729310 9757478 variant2   chr21  9729809 9730866 L1P1        8367  +
chr21  9729310 9757478 variant2   chr21  9730866 9734026 ALR/Alpha   1036  -
chr21  9729310 9757478 variant2   chr21  9734037 9757471 ALR/Alpha   1182  -
chr21  9795588 9796685 variant3   chr21  9795589 9795713 (GAATG)n    308   +
chr21  9795588 9796685 variant3   chr21  9795736 9795894 (GAATG)n    683   +
chr21  9795588 9796685 variant3   chr21  9795911 9796007 (GAATG)n    345   +
chr21  9795588 9796685 variant3   chr21  9796028 9796187 (GAATG)n    756   +
chr21  9795588 9796685 variant3   chr21  9796202 9796615 (GAATG)n    891   +
chr21  9795588 9796685 variant3   chr21  9796637 9796824 (GAATG)n    621   +
$ bedtools groupby -i variantsToRepeats.bed -g 1,2,3 -c 9
chr21 9719758 9729320 6353
chr21 9729310 9757478 14482
chr21 9795588 9796685 3604
$ bedtools groupby -i variantsToRepeats.bed -g 1,2,3 -c 9 -o min
chr21 9719758 9729320 1004
chr21 9729310 9757478 1036
chr21 9795588 9796685 308

参数：

Option	Description
-i	需要进行分类和统计的文件
-g (-grp)	指定哪些列进行分类.指定分类的列可以逗号进行分隔，或者是采用范围的方式(e.g. 1-4)
-c (-opCol)	指定特定的需要进行统计的列. 必须参数
-o (-op)	统计参数的选择，默认是对统计列相同的类别进行加和操作.
-	Valid operations:
-	sum - numeric only
-	count - numeric or text
-	count_distinct - numeric or text
-	min - numeric only
-	max - numeric only
-	mean - numeric only
-	median - numeric only

4.14 intersect

bedtools intersect 命令可以筛选两个基因组特征中的重叠区域。此外，它还设置了一系列参数可以输出不同形式的重复区域，bedtools intersect输入的格式可以是BED/GFF/VCF以及BAM格式。需要注意的是：从版本2.21.0开始，intersect命令的-b参数可以输入多个文件，也就是可以一次性分析一个查询文件（-a参数输入）与多个数据库文件（-b参数输入）的重叠区域。我们现在使用的版本是v2.17.0

示例：

bedtools intersect [OPTIONS] -a <FILE> \
                         -b <FILE1, FILE2, ..., FILEN>

参数：

Option	Description
-a	输入的A文件，格式可以是BAM/BED/GFF/VCF 格式. 每次比较的是A文件与B文件之间的交集
-b	输入的B文件，格式可以是BAM/BED/GFF/VCF格式. 新版本的-b参数可以输入多个文件，寻找多个文件之间的交集
-abam	输入的BAM格式的A文件. 输出BAM文件与 B文件区间的交集. 使用“stdin”代表管道的输入，命令如下: samtools view -b <BAM>/bedtools intersect -abam stdin -b genes.bed
-ubam	输出未压缩的BAM文件。默认输出的压缩格式的BAM文件
-bed	当输入的A文件是BAM格式时(-abam),添加本参数输出的是BED格式. 而不添加本参数，默认输出的是BAM格式
-wa	输出A文件与B文件的重叠区域，并保留A文件原有的输入格式
-wb	输出A文件与B文件的重叠区域，并保留B文件原有的输入格式
-loj	如果A文件与B文件没有重叠区域时，默认不输出不重叠的区域。添加本参数后，不重叠区域也输出，不过以不同形式显示
-wo	添加本参数后，不仅输出A文件与B文件的重叠区域，而且输出重叠区域的大小
-wao	添加本参数后，不仅输出A文件与B文件的重叠区域和重叠区域大小，同时输出不重叠区域，不过以不同形式展示
-u	添加-u 参数只是输出A文件中与B文件有重复区域的区间

intersect 示意图：

intersect.png

4.15 Jaccard

Bedtools intersect 命令可以列举出两个文件中重叠的区间，但是我们经常需要统计的是两个文件中重叠区间的比例来，以此来计算区间的相似程度。Jaccard 命令就统计了重叠区间占整个联合区间的比例。

示例：

bedtools jaccard [OPTIONS] -a <BED/GFF/VCF> -b <BED/GFF/VCF>
$ cat a.bed
chr1  10  20
chr1  30  40

$ cat b.bed
chr1  15   20

$ bedtools jaccard -a a.bed -b b.bed
intersection  union   jaccard n_intersections
5     20      0.25    1

jaccard 示意图：

jaccard.png

4.16 links

该命令会创建一个HTML文件，该文件中的区间链接了UCSC浏览器网站的特定区域，这个工具对于手动查看不同区间的注释信息和其他特征非常有帮助。

示例：

linksBed [OPTIONS] -i <BED/GFF/VCF> > <HTML file>
head genes.bed
chr21 9928613  10012791  uc002yip.1 0  -
chr21 9928613  10012791  uc002yiq.1 0  -
chr21 9928613  10012791  uc002yir.1 0  -
chr21 9928613  10012791  uc010gkv.1 0  -
chr21 9928613  10061300  uc002yis.1 0  -
chr21 10042683 10120796  uc002yit.1 0  -
chr21 10042683 10120808  uc002yiu.1 0  -
chr21 10079666 10120808  uc002yiv.1 0  -
chr21 10080031 10081687  uc002yiw.1 0  -
chr21 10081660 10120796  uc002yix.2 0  -
linksBed -i genes.bed > genes.html

4.17 map

bedtools map命令可以统计B文件中与A文件重叠区域的特征。例如，bedtools map可以计算与A文件有重叠区域的B文件个区间score值的平均值，加和，最小值等指标。

示例：

bedtools map [OPTIONS] -a <bed/gff/vcf> -b <bed/gff/vcf>
$ cat a.bed
chr1        10      20      a1      1       +
chr1        50      60      a2      2       -
chr1        80      90      a3      3       -

$ cat b.bed
chr1        12      14      b1      2       +
chr1        13      15      b2      5       -
chr1        16      18      b3      5       +
chr1        82      85      b4      2       -
chr1        85      87      b5      3       +

$ bedtools map -a a.bed -b b.bed
chr1        10      20      a1      1       +       12
chr1        50      60      a2      2       -       .
chr1        80      90      a3      3       -       5

参数：

Option	Description
-c	统计B文件与A文件重叠区间的特定列，默认是第5列对score值进行统计
-o	对特定列统计的指标：
-	Valid operations:
-	sum - numeric only
-	count - numeric or text
-	count_distinct - numeric or text
-	min - numeric only
-	max - numeric only
-	absmin - numeric only
-	absmax - numeric only
-	mean - numeric only
-	默认是对列的score值进行加和

map 示意图：

map.png

4.18 maskfasta

与getfasta命令相反，bedtools maskfasta 命令可以根据输入的文件的区间遮盖特定的FASTA文件位点，需要注意的是输入的FASTA文件的表头必须和BED文件中染色体的列的名称精确匹配。

示例：

$ bedtools maskfasta [OPTIONS] -fi <input FASTA> -bed <BED/GFF/VCF> -fo <output FASTA>
$ cat test.fa
>chr1
AAAAAAAACCCCCCCCCCCCCGCTACTGGGGGGGGGGGGGGGGGG

$ cat test.bed
chr1 5 10

$ bedtools maskfasta -fi test.fa -bed test.bed -fo test.fa.out

$ cat test.fa.out
>chr1
AAAAANNNNNCCCCCCCCCCGCTACTGGGGGGGGGGGGGGGGGG

参数：

Option	Description
-soft	添加该参数会将遮盖区间内的字母变为小写。不添加该参数，默认将遮盖区间内的字母变为N。
-mc	添加该参数会将遮盖区间内的字母变为输入字符

maskfasta 示意图：

maskfasta.png

4.19 merge

Bedtools merge 命令可以将重叠的区间或者紧邻的区间合并成一个新的区间。

示例：

bedtools merge [OPTIONS] -i <BED/GFF/VCF/BAM>
$ cat A.bed
chr1  100  200
chr1  180  250
chr1  250  500
chr1  501  1000

$ bedtools merge -i A.bed
chr1  100  500
chr1  501  1000

参数：

Option	Description
-s	规定了相同的正负链才能进行合并，默认合并不考虑正负链信息
-S	规定只能选择正链或者负链进行合并‘+’，‘-’。默认合并不考虑正负链信息
-d	规定了两个区间合并的最长距离. 默认的是有重叠区间或者两个区间紧邻或者重叠时才能进行合并

merge 示意图：

merge.png

4.20 multicov

Bedtools mulitcov 命令可以输入多个BAM文件，与BED文件的各区间进行比对，分别计算各个区间重叠的比对数目。

示例：

bedtools multicov [OPTIONS] -bams BAM1 BAM2 BAM3 ... BAMn -bed  <BED/GFF/VCF>
$ cat ivls-of-interest.bed
chr1 0   10000   ivl1
chr1 10000   20000   ivl2
chr1 20000   30000   ivl3
chr1 30000   40000   ivl4

$ bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed ivls-of-interest.bed
chr1 0       10000   ivl1    100 2234    0
chr1 10000   20000   ivl2    123 3245    1000
chr1 20000   30000   ivl3    213 2332    2034
chr1 30000   40000   ivl4    335 7654    0

4.21 window

与bedtools intersect 命令类似，window命令可以寻找A文件和B文件的重叠区域。然而，除了重叠区域，window命令还增加了特定的搜寻范围（默认是增加1000bp），它会将与A文件区间上下游特定范围与B文件区域有重叠的区域也列出。

示例：

$ cat A.bed
chr1  100  200

$ cat B.bed
chr1  500  1000
chr1  1300 2000
$ bedtools window -a A.bed -b B.bed
chr1  100  200  chr1  500  1000
$ cat A.bed
chr1  100  200

$ cat B.bed
chr1  500  1000
chr1  1300 2000
-w 参数列出了增加搜索范围的具体数值
$ bedtools window -a A.bed -b B.bed -w 5000
chr1  100  200  chr1  500   1000
chr1  100  200  chr1  1300  2000

window 示意图：

window.png

4.22 Overlap

Overlap命令可以计算相同行中两个重叠区域的大小

示例：

overlap [OPTIONS] -i <input> -cols s1,e1,s2,e2
windowBed -a A.bed -b B.bed -w 10
chr1  10  20  A  chr1  15  25  B
chr1  10  20  C  chr1  25  35  D
windowBed -a A.bed -b B.bed -w 10 | overlap -i stdin -cols 2,3,6,7
chr1  10  20  A  chr1  15  25  B  5
chr1  10  20  C  chr1  25  35  D  -5

参数：

Option	Description
-i	输入文件.或者用“stdin”代替管道输入.
-cols	指定的列，分别代表两个区间的起始和终止位置。输入的各列必须按指定的方式排列: start1,end1,start2,end2

4.23 random

bedtools random 命令会产生一个BED6格式的随机区间。可以用（-n）参数指定产生区间的数目，可以用（-l）参数指定产生区间的大小。

示例：

bedtools random [OPTIONS] -g <GENOME>
$ bedtools random -g hg19.genome -l 5
chr9  54133731        54133736        1       5       +
chr1  235288830       235288835       2       5       -
chr8  26744718        26744723        3       5       +
chr3  187313616       187313621       4       5       -
chr11 88996846        88996851        5       5       -
chr13 84714855        84714860        6       5       -
chr13 10759738        10759743        7       5       -
chr6  122569739       122569744       8       5       +
chr17 50884025        50884030        9       5       -
chr11 38576901        38576906        10      5       +

参数：

Option	Description
-l	设定产生区间的大小，默认是100bp
-n	设定产生区间的个数.默认是1,000,000个
-seed	设定随机数种子产生特定的区间

random 示意图：

random.png

4.24 reldist

Bedtools reldist 命令可以计算相对距离值，发现一些特定的基因组模式。传统的两组基因组区间相似性比较的方法是基于交叉间隔的数量或者比例。然而，这样的测量方法在很大程度上忽略了两组基因组之间的空间相关性，有些基因区间，尽管间隔一致或者接近，交叉点却很少（例如，增强子和转录起始位点很少重叠，但他们彼此之间的距离比两组随机间隔要近得多）。如果两个集合之间没有空间相关性，那么在相对距离0至0.5之间应该是均匀分布的。然而，如果间隔碰巧比预期近得多，观测到的相对距离的分布就会向较低的相对距离值偏移（如下图）。

reldist1.png

示例：

bedtools reldist [OPTIONS] -a <BED/GFF/VCF> -b <BED/GFF/VCF>
$ bedtools reldist \
   -a data/refseq.chr1.exons.bed.gz \
   -b data/gerp.chr1.bed.gz
reldist  count total fraction
0.00  20629 43422 0.475
0.01  2629  43422 0.061
0.02  1427  43422 0.033
0.03  985 43422 0.023
0.04  897 43422 0.021
0.05  756 43422 0.017
0.06  667 43422 0.015
0.07  557 43422 0.013
0.08  603 43422 0.014

reldist 示意图：

reldist2.png

4.25 shift

Bedtools shift命令可以将区间进行横向移动，相当于
awk '{OFS="\t" print $1,$ 2+<shift>,$3+<shift>}'命令，bedtools shift横移后区间的范围必须限定在染色体的范围之内（不能小于染色体的起始位置，不能大于染色体的终止位置）。

示例：

bedtools shift [OPTIONS] -i <BED/GFF/VCF> -g <GENOME> [-s or (-m and -p)]
$ cat A.bed
chr1 5 100 +
chr1 800 980 -

$ cat my.genome
chr1 1000

$ bedtools shift -i A.bed -g my.genome -s 5
chr1 10 105 +
chr1 805 985 -

$ bedtools shift -i A.bed -g my.genome -p 2 -m -3
chr1 7 102 +
chr1 797 977 -

4.26 shuffle

bedtools shuffle 命令可以随机的变换特征区间在基因组中的位置。还可以提供一个BED/GFF/VCF格式的文件，列出不希望放置的特殊区间。例如，人们不希望一些区间被放置在已知的基因组缺口中。

示例：

bedtools shuffle [OPTIONS] -i <BED/GFF/VCF> -g <GENOME>
$ cat A.bed
chr1  0  100  a1  1  +
chr1  0  1000 a2  2  -

$ cat my.genome
chr1  10000
chr2  8000
chr3  5000
chr4  2000

$ bedtools shuffle -i A.bed -g my.genome
chr4  1498  1598  a1  1  +
chr3  2156  3156  a2  2  -

参数：

Option	Description
-excl	输入的BED文件，指明不希望放置的区间 (e.g., genome gaps).
-incl	输入的BED文件，指明希望放置的区间
-chrom	进行随机变化位置时，指定放置在相同染色体. 默认情况下染色体和区间的位置都是随机选择的
-seed	设定随机数种子产生特定的区间

shuffle 示意图：

shuffle.png

4.27 Slop

与flank功能类似，bedtools slop命令可以根据输入参数扩展输入文件各区间的范围。相当于 awk '{OFS="\t" print $1,$ 2-<slop>,$3+<slop>}'命令，扩展后的区间范围必须限定在染色体的范围之内（不能小于染色体的起始位置，不能大于染色体的终止位置）。

示例：

bedtools slop [OPTIONS] -i <BED/GFF/VCF> -g <GENOME> [-b or (-l and -r)]
$ cat A.bed
chr1 5 100
chr1 800 980

$ cat my.genome
chr1 1000

$ bedtools slop -i A.bed -g my.genome -b 5
chr1 0 105
chr1 795 985

$ bedtools slop -i A.bed -g my.genome -l 2 -r 3
chr1 3 103
chr1 798 983

参数：

Option	Description
-b	规定区间左右两侧分别需要增加的碱基数
-l	规定区间左侧需要增加的碱基数
-r	规定区间右侧需要增加的碱基数

slop 示意图：

slop.png

4.28 sort

bedtools sort 命令可以按照染色体或者其他标准进行排序，默认情况下，bedtools sort 命令首先按照染色体进行升序排列，然后按照起始位置再进行升序排列

示例：

bedtools sort [OPTIONS] -i <BED/GFF/VCF>
cat A.bed
chr1 800 1000
chr1 80  180
chr1 1   10
chr1 750 10000

sortBed -i A.bed
chr1 1   10
chr1 80  180
chr1 750 10000
chr1 800 1000

参数：

Option	Description
-sizeA	按区间大小进行升序排列
-sizeD	按区间大小进行降序排列
-chrThenSizeA	首先根据染色体进行排列,然后按区间大小进行升序排列
-chrThenSizeD	首先根据染色体进行排列,然后按区间大小进行降序排列
-chrThenScoreA	首先根据染色体进行排列,然后按score值进行升序排列
-chrThenScoreD	首先根据染色体进行排列,然后按score值进行降序排列

4.29 Subtract

Bedtools subtract 命令计算区间之间的差集。如果A文件的区间和B文件的区间有重合，那就输出A文件的区间时减去有重合区间的部分。

示例：

bedtools subtract [OPTIONS] -a <BED/GFF/VCF> -b <BED/GFF/VCF>
$ cat A.bed
chr1  10   20
chr1  100  200

$ cat B.bed
chr1  0    30
chr1  180  300

$ bedtools subtract -a A.bed -b B.bed
chr1  100  180

参数：

Option	Description
-f	规定B文件的区间与A文件区间重合的区域占A区间的最小比例，只有重合区间大于这个比例，输出差集时才减去这个重叠区间.默认重叠区间的大小是1bp.
-F	规定B文件的区间与A文件区间重合的区域占B区间的最小比例，只有重合区间大于这个比例，输出差集时才减去这个重叠区间.默认重叠区间的大小是1bp.
-r	规定B文件的区间与A文件区间重合的区域占两个区间的最小比例，只有重合区间大于这个比例，输出差集时才减去这个重叠区间.
-s	规定了取差集必须在相同的正负链之间进行
-S	规定了取差集必须在相反的正负链之间进行
-A	添加此参数，如果A文件的任意区间与B文件的区间有重叠，则不只是重叠的部分被减去，而是整个区间都被减去

subtract 示意图：

subtract.png

4.30 Unionbedg

Bedtools unionbedg 命令可以将多个样本的BED文件组合成一个文件，这样就可以直接来比较多个样本不同区间的基因组特征差异

示例：

bedtools unionbedg [OPTIONS] -i FILE1 FILE2 FILE3 ... FILEn
cat 1.bg
chr1 1000 1500 10
chr1 2000 2100 20

cat 2.bg
chr1 900 1600 60
chr1 1700 2050 50

cat 3.bg
chr1 1980 2070 80
chr1 2090 2100 20

cat sizes.txt
chr1 5000

bedtools unionbedg -i 1.bg 2.bg 3.bg
chr1 900  1000 0  60 0
chr1 1000 1500 10 60 0
chr1 1500 1600 0  60 0
chr1 1700 1980 0  50 0
chr1 1980 2000 0  50 80
chr1 2000 2050 20 50 80
chr1 2050 2070 20 0  80
chr1 2070 2090 20 0  0
chr1 2090 2100 20 0  20

参数：

Option	Description
-header	打印出表头，包括BED文件的染色体，起始位置，终止位置以及定义的BED特征名称
-names	输入一组样品名称列表，添加header参数时名称列表会显示在表头上
-empty	添加基因组文件时（-g），同时输出各样品没有覆盖的区间
-g	添加基因组文件
-filler TEXT	没有覆盖的区域以指定的字符（例如‘N/A’）进行表示. 默认是以0进行表示

5 Bedtools支持的文件格式

5.1 BED格式

完整的BED格式有12列，但是只有前三列是必须的，以下是各列的信息

chrom 染色体的名称
start 从0开始的染色体上的起始位置
end 从1开始的染色体上的终止位置
name 定义的BED的特征名称
score UCSC定的该score值的范围是从0到1000。Bedtools允许任何字符串存储在这个字段中，以便提供更大的注释灵活性。例如该score可以存储科学计数法表示的p值，或者富集平均值
strand 存储正负链信息
thickStart 特征的起始位置（会被bedtools忽略）
thickEnd 特征的终止位置（会被bedtools忽略）
itemRgb RGB值（会被bedtools忽略）
blockCount 间隔或者外显子的数量（会被bedtools忽略）
blockSizes 逗号分隔的间隔大小（会被bedtools忽略）
blockStarts 逗号分隔的间隔的起始位置（会被bedtools忽略）

完整BED格式示例

track name=pairedReads description="Clone Paired Reads" useScore=1
chr22 1000 5000 cloneA 960 + 1000 5000 0 2 567,488, 0,3512
chr22 2000 6000 cloneB 900 - 2000 6000 0 2 433,399, 0,3601

5.2 GFF格式

seqname 序列或者染色体名称例如，“chr1”,“myChrom”
source 特征的来源，包括生成特征的软件，特征来源的数据库
feature 特征的类型例如 "CDS" "start_codon" "stop_codon" and "exon"
start 特征的起始位置
end 特征的终止位置
score 和BED文件的sorce值类似
strand 链的正负值‘+’或‘-’
frame 阅读框，如果是外显子的编码区该值处于0-2之间，如果不是外显子编码区为‘-’
group 所有特征行的分组信息

GFF格式示例

browser position chr22:10000000-10025000
browser hide all
track name=regulatory description="TeleGene(tm) Regulatory Regions" visibility=2
chr22   TeleGene    enhancer    10000000    10001000    500 +   .   touch1
chr22   TeleGene    promoter    10010000    10010100    900 +   .   touch1
chr22   TeleGene    promoter    10020000    10025000    800 -   .   touch2

5.3 BEDPE格式

Bedtools定义了一种新的文件格式（BEDPE）用来处理没有交集的基因组特征文件，如结构变异或者双端序列比对。封闭的BED12格式不能定义染色体间的特性，而且BED12只能表示一条链，对于双端序列比对表示有缺陷，不能很好的展示结构变异。BEDPE文件前10列是必须的，第11列是添加的额外信息。

chrom1 第一个特征的染色体名称
start1 第一个特征的染色体的起始位置
end1 第一个特征的染色体的终止位置
chrom2 第二个特征的染色体的名称
start2 第二个特征的染色体的起始位置
end2 第二个特征的染色体的终止位置
name 定义的BEDPE名称
score UCSC定的该score值的范围是从0到1000。Bedtools允许任何字符串存储在这个字段中，以便提供更大的注释灵活性。例如该score可以存储科学计数法表示的p值，或者富集平均值
strand1 第一个特征的正负链信息‘+’或‘-’
strand2 第二个特征的正负链信息‘+’或‘-’
备选列俩个特征的编辑距离,或者“deletion”,”inversion”等信息

BEDPE格式示例

chr1  10    20    chr5  50    60    a1     30       +    -  0  1
chr9  30    40    chr9  80    90    a2     100      +    -  2  1

5.4 “genome”文件

在使用Bedtools的一些参数（genomecov,complement,slop）时，需要知道BED文件依赖的特定物种的染色体的大小。这时你需要自己构建”gennome”文件，就是简单的罗列出染色体的名称和对应的染色体的大小，并以Tab键进行分隔。

“genome”文件示例

chrI 15072421
chrII 15279323
chrX 17718854
chrM 13794