bedtools 用法大全(一文就够吧)
bedtools等工具号称是可以代替普通的生物信息学数据处理工程师的！我这里用一个专题来讲解它的用法，其实它能实现的需求，我们写脚本都是可以做的，而且我强烈建议正在学编程的小朋友模仿它的各种功能来增强自己的脚本功力。
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）；
好，言归正传，bedtools是非常多的工具的合集，有瑞士军刀的美誉。直接下载二进制版本软件就可以调用全路径来使用，或者把整个文件夹添加到环境变量。
首发于 生信技能树：http://www.biotrainee.com/thread-153-1-3.html
第一个功能 genomecov
我们先看第一个功能，把alignment的结果文件转为bedgraph格式文件。 不过这个功能用处不是很大。
参考：http://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html
要实现这个功能，需要用bedtools的genomecov小命令， 有两种方法可以调用：
bedtools genomecov [OPTIONS] [-i|-ibam] -g (iff. -i)
genomeCoverageBed [OPTIONS] [-i|-ibam] -g (iff. -i)
这个命令本身并不是设计来做格式转换的，bam2bedgraph也只是其中的一个小功能而已，需要加上-bg参数，就可以Report depth in BedGraph format. For details, see: http://genome.ucsc.edu/goldenPath/help/bedgraph.html
大家观摩我下面给出的测试例子，就明白该功能如何使用了
bedtools genomecov -bg -i E001-H3K4me1.tagAlign -g mygenome.txt >E001-H3K4me1.bedGraph
bedtools genomecov -bg -i E001-Input.tagAlign -g mygenome.txt >E001-Input.bedGraph
nohup bedtools genomecov -bg -ibam BAF180_CT10.unique.sorted.bam >BAF180_CT10.bedGraph &
nohup bedtools genomecov -bg -ibam BAF180_CT22AM.unique.sorted.bam >BAF180_CT22AM.bedGraph &
nohup bedtools genomecov -bg -ibam BAF180_CT22.unique.sorted.bam >BAF180_CT22.bedGraph &
nohup bedtools genomecov -bg -ibam inputCT10sonication.unique.sorted.bam >inputCT10sonication.bedGraph &
首先alignment的文件必须是sort的，然后如果是bed格式的比对文件，用-i 参数来指定输入文件，需要加入参考基因组的染色体大小记录文件(mygenome.txt )，如果是bam格式的比对文件，用-ibam指定输入文件，而且不需要参考基因组的染色体大小记录文件。
第二个功能 multicov
然后看看第二个功能,对RNA-seq的比对文件中的比对到各个基因的reads进行计数。
参考: http://www.bio-info-trainee.com/745.html
实现这个功能，用的bedtools的multicov 这个小命令
# 例子：
bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed ivls-of-interest.bed
# ivls-of-interest.bed这个文件是必须的，可能需要自己制作，其实用gtf文件也可以的，如下：
chr1 0 10000 ivl1
chr1 10000 20000 ivl2
chr1 20000 30000 ivl3
chr1 30000 40000 ivl4
输出结果前三列是坐标，第四列是基因名，跟我们的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
可以看到，它实现的需求，跟htseq这个软件差不多。各种区别，大家可以自己取探索。
第三个功能 getfasta
接着第三个功能，根据坐标区域来从基因组里面提取fasta序列
参考：# BED/GFF/VCF +reference --> fasta
bedtools getfasta -fi ~/biosoft/bowtie/hg19_index/hg19.fa -bed ../macs14_results/highQuality_summits.bed -fo highQuality.fa
bedtools getfasta -fi ~/biosoft/bowtie/hg19_index/hg19.fa -bed ../macs14_results/highQuality_peaks.bed -fo highQuality.fa
我的例子脚本里面用的是bed格式来记录坐标区域，参考基因组用-fi参数指定具体位置，输出的fasta序列文件用-fo参数指定
第四个功能区域注释 intersect
首发于生信技能树论坛：http://www.biotrainee.com/thread-1700-1-2.html
下载的CNV都是基于基因组区域的，比如1号染色体的61735起始坐标到1510801终止坐标。在IGV里面倒是可以看出一些pattern，但是人们感兴趣的往往是这些位置上面到底有哪些基因。接下来就可以对基因进行各种下游分析。
既然是对基因组片段做基因注释，那么首先就需要拿到基因的坐标信息咯，我是在gencode数据库里面下载，然后解析成下面的bed格式的，如下：
head ~/reference/gtf/gencode/protein_coding.hg19.position
chr1 69091 70008 OR4F5
chr1 367640 368634 OR4F29
chr1 621096 622034 OR4F16
chr1 859308 879961 SAMD11
chr1 879584 894689 NOC2L
chr1 895967 901095 KLHL17
chr1 901877 911245 PLEKHN1
chr1 910584 917473 PERM1
chr1 934342 935552 HES4
chr1 936518 949921 ISG15
然后要把我们下载的CNV文本文件，转为bed格式的，就是把列的顺序调换一下：
head Features.bed
chr1 3218610 95674710 TCGA-3C-AAAU-10A-01D-A41E-01 53225 0.0055
chr1 95676511 95676518 TCGA-3C-AAAU-10A-01D-A41E-01 2 -1.6636
chr1 95680124 167057183 TCGA-3C-AAAU-10A-01D-A41E-01 24886 0.0053
chr1 167057495 167059336 TCGA-3C-AAAU-10A-01D-A41E-01 3 -1.0999
chr1 167059760 181602002 TCGA-3C-AAAU-10A-01D-A41E-01 9213 -8e-04
chr1 181603120 181609567 TCGA-3C-AAAU-10A-01D-A41E-01 6 -1.2009
chr1 181610685 201473647 TCGA-3C-AAAU-10A-01D-A41E-01 12002 0.0055
chr1 201474400 201474544 TCGA-3C-AAAU-10A-01D-A41E-01 2 -1.4235
chr1 201475220 247813706 TCGA-3C-AAAU-10A-01D-A41E-01 29781 -4e-04
命令很简单，如下：
bedtools intersect -a Features.bed -b ~/reference/gtf/gencode/protein_coding.hg19.position \
-wa -wb | bedtools groupby -i - -g 1-4 -c 10 -o collapse
注释结果如下：可以看到，每个CNV片段都注释到了对应的基因，有些特别大的片段，会被注释到非常多的基因。
chr8 42584924 42783715 TCGA-5T-A9QA-01A-11D-A41E-01 CHRNB3,CHRNA6,THAP1,RNF170,HOOK3
chr8 42789728 42793594 TCGA-5T-A9QA-01A-11D-A41E-01 HOOK3
chr8 42797957 42933372 TCGA-5T-A9QA-01A-11D-A41E-01 RP11-598P20.5,FNTA,HOOK3
chr8 70952673 70964372 TCGA-5T-A9QA-01A-11D-A41E-01 PRDM14
chr10 42947970 43833200 TCGA-5T-A9QA-01A-11D-A41E-01 BMS1,RET,RASGEF1A,ZNF33B,CSGALNACT2
chr10 106384615 106473355 TCGA-5T-A9QA-01A-11D-A41E-01 SORCS3
chr10 106478366 107298256 TCGA-5T-A9QA-01A-11D-A41E-01 SORCS3
chr10 117457285 117457859 TCGA-5T-A9QA-01A-11D-A41E-01 ATRNL1
chr11 68990173 69277078 TCGA-5T-A9QA-01A-11D-A41E-01 MYEOV
chr11 76378708 76926535 TCGA-5T-A9QA-01A-11D-A41E-01 LRRC32,B3GNT6,OMP,TSKU,MYO7A,ACER3,CAPN5
最常用的 intersect 的8个案例
用来求两个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
4.bedtools merge
用于合并位于同一个bed/gff/vcf 文件中的重叠区域。
Bedtools merge [OPTION] –i
-s 必须相同(正负)链的区域才合并（软件默认不考虑正负链特征）
-n 报告合并的区域数量，没有被合并则1
-d 两个独立区域间距小于（等于）该值时将被合并为一个区域
-nms 报告被合并区域的名称
-scores 报告几个被合并特征区域的scores
其他小功能
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).
有很多组织单位给bedtools写说明书：
http://bedtools.readthedocs.io/en/latest/index.html
http://quinlanlab.org/tutorials/bedtools/bedtools.html
猜你喜欢
生信基础知识100讲
还有更多文章，请移步公众号阅读
▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。
▼ 如果你是初学者，请关注下面的生信菜鸟团，了解生信基础名词，概念，扎实的打好基础，争取早日入门。