bedtools用法大全一文就够吧

1、bedtools 用法大全 (一文就够吧 )bedtools 等工具号称是可以代替普通的生物信息学数据 处理工程师的！我这里用一个专题来讲解它的用法，其实它 能实现的需求，我们写脚本都是可以做的，而且我强烈建议 正在学编程的小朋友模仿它的各种功能来增强自己的脚本 功力。 BEDTools 是可用于 genomic features 的比较，相关 操作及进行注释的工具。而 genomic features 通常使用 Browser Extensible Data (BED) 或者 General Feature Format (GFF) 文件表示，用 UCSC Genome Browser 进行

2、 可视化比较。 bedtools 总共有二三十个工具 /命令来处理基因 组数据。比较典型而且常用的功能举例如下： 格式转换， bam 转 bed ( bamToBed )， bed 转其他格式( bedToBam ， bedToIgv )； 对基因组坐标的逻辑运算，包括：交集(intersectBed , windowBed ), ” 邻集“(osestBed ), 补集( complementBed )，并集( mergeBed )，差集( subtractBed );计算覆盖度( coverage )( coverageBed ， genomeCoverageBed )；好，言归正传， b

3、edtools 是非常多的工具的合集，有瑞士军 刀的美誉。直接下载二进制版本软件就可以调用全路径来使用，或者把整个文件夹添加到环境变量。首发于生信技能第一个树： 功能 genomecov 我们先看第一个功能，把 alignment 的结 果文件转为 bedgraph 格式文件。 不过这个功能用处不是很http:/bedtools.readthedocs.io/en/latest/content/tools/geno mecov.html 要实现这个功能，需要用 bedtools 的genomecov 小命令， 有两种方法可以调用： bedtools genomecov OPTIONS -i|-

4、ibam -g (iff. -i)genomeCoverageBed OPTIONS -i|-ibam -g (iff. -i)这个命令本身并不是设计来做格式转换的，bam2bedgraph也只是其中的一个小功能而已，需要加上 -bg 参数，就可以Report depth in BedGraph format. For details, see: /goldenPath/help/bedgraph.html家观摩我下面给出的测试例子，就明白该功能如何使用了bedtools genomecov -bg -i E001-H3K4me1.tagAlign -

5、g mygenome.txt E001-H3K4me1.bedGraph bedtools genomecov -bg -i E001-Input.tagAlign -g mygenome.txt E001-Input.bedGraph nohup bedtools genomecov -bg -ibamBAF180_CT10.unique.sorted.bam BAF180_CT10.bedGr aph &nohup bedtools genomecov-bg -ibamBAF180_CT22AM.unique.sorted.bam BAF180_CT22AM. bedGraph &nohu

6、p bedtools genomecov-bg -ibamBAF180_CT22.unique.sorted.bam BAF180_CT22.bedGr aph &nohup bedtools genomecov-bg -ibaminputCT10sonication.unique.sorted.bam inputCT10sonic ation.bedGraph &首先 alignment 的文件必须是 sort 的，然后如果是 bed 格式 的比对文件，用 -i 参数来指定输入文件，需要加入参考基因 组的染色体大小记录文件 (mygenome.txt ) ，如果是 bam 格 式的比对文件，

7、用 -ibam 指定输入文件，而且不需要参考基 因组的染色体大小记录文件。 第二个功能 multicov 然后看看 第二个功能 ,对 RNA-seq 的比对文件中的比对到各个基因的 reads 进行计数。参考 : http:/www.bio-info- 实现这个功能， 用 的 bedtools 的 multicov 这个小命令 # 例子： bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed ivls-of-interest.bed# ivls-of-interest.bed 这个文件是必须的， 可能需要自己制作， 其实用 gtf 文件

8、也可以的，如下：chr1 0 10000 ivl1chr1 1000020000ivl2chr1 2000030000ivl3chr1 3000040000ivl4输出结果前三列是坐标，第四列是基因名，跟我们的 bed 文 件一样，只是最后三列是三个样本的计数，是添加上来的！chr1010000ivl1100 22340chr11000020000ivl2123 32451000chr12000030000ivl3213 23322034chvl4335 76540可以看到，它实现的需求，跟 htseq 这个软件差不多。各种 区别， 大家可以自己取探索。 第三个功能

9、getfasta 接着第三 个功能，根据坐标区域来从基因组里面提取 fasta 序列参考： # BED/GFF/VCF +reference - fastabedtools getfasta -fi /biosoft/bowtie/hg19_index/hg19.fa -bed ./macs14_results/highQuality_summits.bed-fohighQuality.fa bedtools getfasta -fi /biosoft/bowtie/hg19_index/hg19.fa -bed ./macs14_results/highQuality_peaks.bed-f

10、o highQuality.fa我的例子脚本里面用的是 bed 格式来记录坐标区域， 参考基 因组用 -fi 参数指定具体位置， 输出的 fasta 序列文件用 -fo 参 数指定第四个功能区域注释 intersect 首发于生信技能树论 坛： 下载的 CNV 都是基于基因组区域的，比如 1 号染色体的 61735 起 始坐标到 1510801 终止坐标。在 IGV 里面倒是可以看出一 些 pattern ，但是人们感兴趣的往往是这些位置上面到底有哪 些基因。接下来就可以对基因进行各种下游分析。既然是对 基因组片段做基因注释，那么首先就需要拿到基因的坐标信 息咯，我是在 gencode 数据库

11、里面下载， 然后解析成下面的 bed 格式的，如下： head/reference/gtf/gencode/protein_coding.hg19.positionchr16909170008OR4F5chr1367640368634OR4F29chr1621096622034OR4F16chr1859308879961SAMD11chr1879584894689NOC2Lchr1895967901095KLHL17chr1901877911245PLEKHN1chr1910584917473PERM1chr1934342935552HES4chr1936518949921ISG15然后要把我

12、们下载的CNV文本文件，转为 bed 格式的，就是把列的顺序调换一下： head Features.bedchr13218610 95674710TCGA-3C-AAAU-10A-01D-A41E-01 chr19567651195676518TCGA-3C-AAAU-10A-01D-A41E-01chr195680124167057183TCGA-3C-AAAU-10A-01D-A41E-01chr1167057495167059336TCGA-3C-AAAU-10A-01D-A41E-01chr1167059760181602002TCGA-3C-AAAU-10A-01D-A41E-01c

13、hr1181603120181609567TCGA-3C-AAAU-10A-01D-A41E-01chr1181610685201473647TCGA-3C-AAAU-10A-01D-A41E-01chr1201474400201474544TCGA-3C-AAAU-10A-01D-A41E-01chr1201475220247813706TCGA-3C-AAAU-10A-01D-A41E-0153225 0.00552 -1.663624886 0.00533 -1.09999213 -8e-046 -1.200912002 0.00552 -1.423529781 -4e-04命令很简单，

14、如下： 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 42783715TCGA-5T-A9QA-01A-11D-A41E-01CHRNB3,CHRNA6,THAP1,RNF170,HOOK3chr8427897284

15、2793594TCGA-5T-A9QA-01A-11D-A41E-01 HOOK3 chr84279795742933372TCGA-5T-A9QA-01A-11D-A41E-01RP11-598P20.5,FNTA,HOOK3chr8 70952673 70964372TCGA-5T-A9QA-01A-11D-A41E-01 PRDM14 chr10 42947970 43833200TCGA-5T-A9QA-01A-11D-A41E-01BMS1,RET,RASGEF1A,ZNF33B,CSGALNACT2chr10106384615106473355TCGA-5T-A9QA-01A-11

16、D-A41E-01SORCS3chr10106478366107298256TCGA-5T-A9QA-01A-11D-A41E-01SORCS3chr10117457285117457859TCGA-5T-A9QA-01A-11D-A41E-01ATRNL1chr116899017369277078TCGA-5T-A9QA-01A-11D-A41E-01MYEOVchr117637870876926535TCGA-5T-A9QA-01A-11D-A41E-01 LRRC32,B3GNT6,OMP,TSKU,MYO7A,ACER3,CAPN5 最常用的 intersect 的 8 个案例用来求两

17、个 BED 或者 BAM 文件中的 overlap ， overlap 可以进行自定义是整个 genome features 的 overlap 还是局部。 加 -wa 参数可以报告出原始 的在 A 文件中的 feature ，加 -wb 参数可以报告出原始的在 B 文件中的 feature, 加 -c 参数可以报告出两个文件中的 overlap 的 feature 的数量，参数 -s 可以得到忽略 strand 的 overlap ，具体案例如下：案例一：包含着染色体位置的两个 文件， 分别记为 A 文件和 B 文件。 分别来自于不同文件的染 色体位置的交集是什么？$ cat A.bed c

18、hr1 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.bedchr1 10 20 chr1 30 40$ cat B.bedchr1 15 25$ bedtools intersect -a A.bed -b B.bed -wachr1 10 20 案例三：包含着染色体位置的两个文件，分别记为A文件和B

19、文件。求A文件中染色体位置与文件B中染色体位置的交集，以及对应的文件 B 中的染色体位置 .$ cat A.bedchr1 10 20chr1 30 40$ cat B.bedchr1 15 25$ bedtools intersect -a A.bed -b B.bed -wbchr1 15 20 chr1 15 25 案例四（经用） ： 包含着染色体位置 的两个文件，分别记为A 文件和 B 文件。求对于 A 文件的染色体位置是否与文件B 中的染色体位置有交集。如果有交集，分别输入 A 文件的染色体位置和 B 文件的染色体位置； 如果没有交集， 输入 A 文件的染色体位置并以 . -1 -1

20、 补齐文 件。$ cat A.bed chr1 10 20chr1 30 40 $ cat B.bed chr1 15 25$ bedtools intersect -a A.bed -b B.bed -lojchr1 10 20 chr1 15 25chr1 30 40 . -1 -1 案例五： 包含着染色体位置的两个文件， 分别记为 A 文件和 B 文件。对于 A 文件中染色体位置，如 果和 B 文件中染色体位置有 overlap, 则输出在 A 文件中染色 体位置和在 B 文件中染色体位置，以及 overlap 的长度 . $ cat A.bed chr1 10 20chr1 30 40

21、 $ cat B.bed chr1 15 20 chr1 18 25 $ bedtools intersect -a A.bed -b B.bed -wo chr1 10 20 chr1 15 20 5chr1 10 20 chr1 18 25 2 案例六： 包含着染色体位置的两个 文件，分别记为 A 文件和 B 文件。对于 A 文件中染色体位 置，如果和 B 文件中染色体位置有 overlap, 则输出在 A 文件 中染色体位置和在 B 文件中染色体位置，以及 overlap 的长 度；如果和 B 文件中染色体位置都没有 overlap, 则用 . -1-1补齐文件$ cat A.bed c

22、hr1 10 20chr1 30 40$ cat B.bedchr1 15 20chr1 18 25$ bedtools intersect -a A.bed -b B.bed -waochr1 10 20 chr1 15 20 5 chr1 10 20 chr1 18 25 2chr1 30 40 . -1 -1 案例七： 包含着染色体位置的两个文件， 分别记为 A 文件和 B 文件。对于 A 文件中染色体位置，输 出在A文件中染色体位置和有多少B文件染色体位置与之有overlap.$ cat A.bedchr1 10 20 chr1 30 40$ cat B.bedchr1 15 20ch

23、r1 18 25$ bedtools intersect -a A.bed -b B.bed -c chr1 10 20 2chr1 30 40 0 案例八 （常用）：包含着染色体位置的两个文件， 分别记为 A 文件和 B 文件。对于 A 文件中染色体位置，输 出在 A 文件中染色体位置和与 B 文件染色体位置至少有 X% 的 overlap 的记录。$ cat A.bedchr1 100 200$ cat B.bedchr1 130 201chr1 180 220$ bedtools intersect -a A.bed -b B.bed -f 0.50 -wa -wbchr1 100 200 chr1 130 2014.bedtools merge 用于合并位于同一个 bed/gff/vcf 文件中的重叠区域。Bedtools merge OPTION -s 必须相同 （正负 ）链的区域才合并（软件默认不考虑正负链 特征）-n 报告合并的区域数量，没有被合并则 1-d 两个独立区域间距小于 （等于） 该值时将被合并为一个区域-nms 报告被合并区域的名称-scores 报告几个被合并