bedtools开发的目的是为了快速，灵活的比较大量的基因组特征（genomic features）。而genomic features通常使用Browser Extensible Data (BED) 或者 General Feature Format (GFF)文件表示，用UCSC Genome Browser进行可视化比较。

例如,bedtools可以进行取intersect（交集）, merge（并集）, count（计数）, complement（补集），以及用来对广泛使用的基因组文件格式，例如BAM, BED, GFF/GTF, VCF等进行基因组区间的转换。单个的工具设计的目的是应对简单的任务，复杂的分析能通过组合多个bedtools工具操作实现。同时，该工具允许控制输出结果的呈现形式。最初的bedtools版本支持单独的6列BED文件。但是，如今增加了对序列比对BAM文件的支持。以及GFF文件的特征，BED文件。以及VCF文件。这些工具是相当快速的，并且即使是大的数据集也可以在数秒内完成任务。

我们用bedtools都可以做些啥?

bedtools总共有二三十个工具/命令来处理基因组数据。比较典型而且常用的功能举例如下：格式转换，bam转bed（bamToBed），bed转其他格式（bedToBam，bedToIgv）；对基因组坐标的逻辑运算，包括：交集（intersectBed，windowBed），”邻集“（closestBed），补集（complementBed），并集（mergeBed），差集（subtractBed）;计算覆盖度（coverage）（coverageBed，genomeCoverageBed）；
此外，还有一些强大而实用的工具（shuffleBed，groupBy，annotateBed，……）

Utility	Description
annotate	Annotate coverage of features from multiple files.
bamtobed	Convert BAM alignments to BED (& other) formats.
bamtofastq	Convert BAM records to FASTQ records.
bed12tobed6	Breaks BED12 intervals into discrete BED6 intervals.
bedpetobam	Convert BEDPE intervals to BAM records.
bedtobam	Convert intervals to BAM records.**
closest	Find the closest, potentially non-overlapping interval.
cluster	Cluster (but don’t merge) overlapping/nearby intervals.
complement	Extract intervals not represented by an interval file.
coverage	Compute the coverage over defined intervals.
expand	Replicate lines based on lists of values in columns.
flank	Create new intervals from the flanks of existing intervals.
genomecov	Compute the coverage over an entire genome.
getfasta	Use intervals to extract sequences from a FASTA file.
groupby	Group by common cols. & summarize oth. cols. (~ SQL “groupBy”)
igv	Create an IGV snapshot batch script.
intersect	Find overlapping intervals in various ways.
jaccard	Calculate the Jaccard statistic b/w two sets of intervals.
links	Create a HTML page of links to UCSC locations.
makewindows	Make interval “windows” across a genome.
map	Apply a function to a column for each overlapping interval.
maskfasta	Use intervals to mask sequences from a FASTA file.
merge	Combine overlapping/nearby intervals into a single interval.
multicov	Counts coverage from multiple BAMs at specific intervals.
multiinter	Identifies common intervals among multiple interval files.
nuc	Profile the nucleotide content of intervals in a FASTA file.
overlap	Computes the amount of overlap from two intervals.
pairtobed	Find pairs that overlap intervals in various ways.
pairtopair	Find pairs that overlap other pairs in various ways.
random	Generate random intervals in a genome.
reldist	Calculate the distribution of relative distances b/w two files.
shift	Adjust the position of intervals.
shuffle	Randomly redistribute intervals in a genome.
slop	Adjust the size of intervals.
sort	Order the intervals in a file.
subtract	Remove intervals based on overlaps b/w two files.
tag	Tag BAM alignments based on overlaps with interval files.
unionbedg	Combines coverage intervals from multiple BEDGRAPH files.
window	Find overlapping intervals within a window around an interval.

BEDTools suite使用详细

bedtools官网：

http://bedtools.readthedocs.io/en/latest/

bedtools使用说明：

http://quinlanlab.org/tutorials/bedtools/bedtools.html#bedtools-merge

BEDTools主要使用BED格式的前三列,即：

chrom: 染色体信息
start: genome feature的起始位点，从0开始
end: genome feature的终止位点，至少为1

一般常用物种的genome file在BEDTools安装目录的/genome里面

BEDPE格式是其自定义的一种新的格式，为了简洁的描述不连续的genome features，例如结构变异和双端测序比对

注意：

start1和start2起始坐标第一个碱基都为0，所以start=9, end=20表示碱基跨度是从第10位到第20位
chrom1或者chrom2用.表示unknown;start1，end1,start2,end2用-1表示unknown

(1)intersect

可以计算两个或者多个BED/BAM/VCF/GFF文件中基因组坐标位置的交集(overlap)，根据参数不同，可以得到不同的结果。

两个BED文件比较图示

一对多比较图示

语法：

bedtools intersect  -a <bed/gff/vcf/bam> -b <bed/gff/vcf/bam> [OPTIONS]

-wa参数可以报告出原始的在A文件中的feature
-wb参数可以报告出原始的在B文件中的feature
-c参数可以报告出两个文件中的overlap的feature的数量
-wo 返回overlap碱基数
-v 返回非overlap区间
-s 相同链上的feature

当用bedtools intersect 处理大文件时比较耗内存，有效的方法是对A和B文件按照染色体名字(chromosome)和位置(position)排序(sort -k1,1 -k2,2n),然后用-sorted参数重新intersect

案例

注意，自己生成测试bed文件，都必须用tab键分割，否则会报错！！

案例一：包含着染色体位置的两个文件，分别记为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)merge

用于合并位于同一个bed/gff/vcf 文件中有overlap或者距离在一定范围内的相邻区间，距离可由参数(-d)定义。需要注意的是，做合并之前需要先对bed文件做根据染色体排序，可以用bedtoolssort命令实现。

图示

语法

bedtools merge [OPTIONS] -i <bed/gff/vcf>

案例

$ cat A.bed

chr2 10 20

chr1 30 40

chr1 15 20

chr1 18 25

排序：

sort -k1,1 -k2,2n A.bed > A.sort.bed

$ cat A.sort.bed 

chr1 15 20

chr1 18 25

chr1 30 40

chr2 10 20

案例一：取并集

bedtools merge -i A.sort.bed

chr1 15 25

chr1 30 40

chr2 10 20
案例二：计算重叠区间的个数,-i 指定统计的列，-o指定操作5

bedtools merge -i exons.bed -c 1 -o count

chr1 15 25 2

chr1 30 40 1

chr2 10 20 1
案例三：-d 两个独立区域间距小于（等于）该值时将被合并为一个区域；-o collapse显示合并了哪些标签

$ bedtools merge -i A.sort.bed -d 5 -c 1 -o count，collapse

chr1 15 40 3 chr1,chr1,chr1

chr2 10 20 1 chr2

(3)complement：返回基因组非覆盖区（用途，比如多轮设计panel)

图示

语法

bedtools complement -i <BED/GFF/VCF> -g <genome files>

(4)genomecov：染色体和全基因组覆盖度计算

要求：单个输入bed文件（-i指定）和genome files；如果输入为bam(-ibam指定)文件，则不需要genome files

图示

语法

bedtools genomecov [OPTIONS] -i <bed/gff/vcf> -g <genome>

案例

$ cat ranges-cov-sorted.bed
chr1    4       9
chr1    1       6
chr1    8       19
chr1    25      30
chr2    0       20

$ cat cov.txt  （染色体及每条染色体总碱基数）
chr1    30
chr2    20

bedtools genomecov -i ranges-cov-sorted.bed -g cov.txt
chr1    0       7       30      0.233333 1
chr1    1       20      30      0.666667
chr1    2       3       30      0.1
chr2    1       20      20      1 2
genome  0       7       50      0.14 3
genome  1       40      50      0.8
genome  2       3       50      0.06
#name 覆盖次数 覆盖碱基数 总碱基数   覆盖度
#同时计算单染色体和全基因组覆盖度

ranges-cov.bed文件需提前排序sort -k1,1 ranges-cov.bed > ranges-cov-sorted.bed
-bg参数可得到每个碱基的覆盖度。

(5)coverage 计算染色体给定区间覆盖度，输入可以是 BAM 文件

$ 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
#name 覆盖次数 覆盖碱基数 总碱基数   覆盖度

(6)getfasta：提取序列

提取指定位置的 DNA 序列，也是很好用的一个功能，反向互补链也可以提，不用自己写脚本提了


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

注意取的是（5,10]，区间左开右闭

语法

bedtools getfasta [OPTIONS] -fi <fasta> -bed <bed/gff/vcf> -fo <fasta>

案例

要求：基因组fasta文件（-fi指定）,提取区间BED/GTF/GFF/VCF文件(-bed指定),输出文件FASTA（-fo 指定）

bedtools getfasta -fi Mus_musculus.GRCm38.75.dna_rm.toplevel_chr1.fa -bed mm_GRCm38_3kb_promoters.gtf -fo mm_GRCm38_3kb_promoters.fasta

扩展：

提取序列之samtools（速度较快）

#首先建立fai索引文件（第一列为染色体名字，第二列为序列碱基数）

samtools faidx Mus_musculus.GRCm38.75.dna.chromosome.8.fa

#序列提取，多提取区间空格隔开

samtools faidx Mus_musculus.GRCm38.75.dna.chromosome.8.fa \
 8:123407082-123410744 8:123518835-123536649

>8:123407082-123410744
GAGAAAAGCTCCCTTCTTCTCCAGAGTCCCGTCTACCCTGGCTTGGCGAGGGAAAGGAAC
CAGACATATATCAGAGGCAAGTAACCAAGAAGTCTGGAGGTGTTGAGTTTAGGCATGTCT
[...]

>8:123518835-123536649
TCTCGCGAGGATTTGAGAACCAGCACGGGATCTAGTCGGAGTTGCCAGGAGACCGCGCAG
CCTCCTCTGACCAGCGCCCATCCCGGATTAGTGGAAGTGCTGGACTGCTGGCACCATGGT
[...]

(7)nuc: 计算GC含量即各碱基数

语法

bedtools nuc [OPTIONS] -fi <fasta> -bed <bed/gff/vcf>

Options:

-fi     输入FASTA文件

-bed     提取区间BED/GTF/GFF/VCF文件(-bed指定)

案例

bedtools nuc -fi hg19.fa -bed CDS.bed

输出结果解释：在原bed文件每行结尾增加以下几列

Output format:
The following information will be reported after each BED entry:

1) %AT content
2) %GC content
3) Number of As observed
4) Number of Cs observed
5) Number of Gs observed
6) Number of Ts observed
7) Number of Ns observed
8) Number of other bases observed
9) The length of the explored sequence/interval.
10) The seq. extracted from the FASTA file. (opt., if -seq is used)
11) The number of times a user's pattern was observed.
    (opt., if -pattern is used.)