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
Bedtools annotate 命令用于统计一个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 示意图:
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 示意图：
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 示意图：
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 示意图：
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 示意图：
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 示意图：
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 示意图：
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 示意图：
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 示意图：
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 示意图：
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 示意图：
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 示意图：
4.24 reldist
Bedtools reldist 命令可以计算相对距离值，发现一些特定的基因组模式。传统的两组基因组区间相似性比较的方法是基于交叉间隔的数量或者比例。然而，这样的测量方法在很大程度上忽略了两组基因组之间的空间相关性，有些基因区间，尽管间隔一致或者接近，交叉点却很少（例如，增强子和转录起始位点很少重叠，但他们彼此之间的距离比两组随机间隔要近得多）。如果两个集合之间没有空间相关性，那么在相对距离0至0.5之间应该是均匀分布的。然而，如果间隔碰巧比预期近得多，观测到的相对距离的分布就会向较低的相对距离值偏移（如下图）。
示例：
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 示意图：
4.25 shift
Bedtools shift命令可以将区间进行横向移动，相当于
awk '{OFS="\t" print 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 示意图：
4.27 Slop
与flank功能类似，bedtools slop命令可以根据输入参数扩展输入文件各区间的范围。相当于 awk '{OFS="\t" print 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 示意图：
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 示意图：
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
