bedtools用法大全一文就够吧.docx

1、bedtools用法大全一文就够吧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），”邻集“（closestBed），补集（complementBed），并集（mergeBed），差集（subtractBed）;计算覆盖度（coverage）（coverageBed，genomeCoverageBed）；好，言归正传，bedtools是非常多的工具的合集，有瑞士军刀的美誉。直接下载二进制版本软

3、件就可以调用全路径来使用，或者把整个文件夹添加到环境变量。首发于 生信技能树： 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

4、 -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.bedGraphbedtools genomec

5、ov -bg -i E001-Input.tagAlign -g mygenome.txt E001-Input.bedGraphnohup 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.sor

6、ted.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-se

7、q的比对文件中的比对到各个基因的reads进行计数。参考: http:/www.bio-info- 这个小命令# 例子：bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed ivls-of-interest.bed# ivls-of-interest.bed这个文件是必须的，可能需要自己制作，其实用gtf文件也可以的，如下：chr1 0 10000 ivl1chr1 10000 20000 ivl2chr1 20000 30000 ivl3chr1 30000 40000 ivl4输出结果前三列是坐标，第四列是基因名，跟我们的bed文件

8、一样，只是最后三列是三个样本的计数，是添加上来的！chr1 0 10000 ivl1 100 2234 0chr1 10000 20000 ivl2 123 3245 1000chr1 20000 30000 ivl3 213 2332 2034chr1 30000 40000 ivl4 335 7654 0可以看到，它实现的需求，跟htseq这个软件差不多。各种区别，大家可以自己取探索。第三个功能 getfasta接着第三个功能，根据坐标区域来从基因组里面提取fasta序列参考：# BED/GFF/VCF +reference - fastabedtools getfasta -fi /bi

9、osoft/bowtie/hg19_index/hg19.fa -bed ./macs14_results/highQuality_summits.bed -fo highQuality.fabedtools getfasta -fi /biosoft/bowtie/hg19_index/hg19.fa -bed ./macs14_results/highQuality_peaks.bed -fo highQuality.fa我的例子脚本里面用的是bed格式来记录坐标区域，参考基因组用-fi参数指定具体位置，输出的fasta序列文件用-fo参数指定第四个功能区域注释 intersect首发于生

10、信技能树论坛： /reference/gtf/gencode/protein_coding.hg19.positionchr1 69091 70008 OR4F5chr1 367640 368634 OR4F29chr1 621096 622034 OR4F16chr1 859308 879961 SAMD11chr1 879584 894689 NOC2Lchr1 895967 901095 KLHL17chr1 901877 911245 PLEKHN1chr1 910584 917473 PERM1chr1 934342 935552 HES4chr1 936518 949921 ISG

11、15然后要把我们下载的CNV文本文件，转为bed格式的，就是把列的顺序调换一下：head Features.bed chr1 3218610 95674710 TCGA-3C-AAAU-10A-01D-A41E-01 53225 0.0055chr1 95676511 95676518 TCGA-3C-AAAU-10A-01D-A41E-01 2 -1.6636chr1 95680124 167057183 TCGA-3C-AAAU-10A-01D-A41E-01 24886 0.0053chr1 167057495 167059336 TCGA-3C-AAAU-10A-01D-A41E-01

12、 3 -1.0999chr1 167059760 181602002 TCGA-3C-AAAU-10A-01D-A41E-01 9213 -8e-04chr1 181603120 181609567 TCGA-3C-AAAU-10A-01D-A41E-01 6 -1.2009chr1 181610685 201473647 TCGA-3C-AAAU-10A-01D-A41E-01 12002 0.0055chr1 201474400 201474544 TCGA-3C-AAAU-10A-01D-A41E-01 2 -1.4235chr1 201475220 247813706 TCGA-3C-

13、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,CH

14、RNA6,THAP1,RNF170,HOOK3chr8 42789728 42793594 TCGA-5T-A9QA-01A-11D-A41E-01 HOOK3chr8 42797957 42933372 TCGA-5T-A9QA-01A-11D-A41E-01 RP11-598P20.5,FNTA,HOOK3chr8 70952673 70964372 TCGA-5T-A9QA-01A-11D-A41E-01 PRDM14chr10 42947970 43833200 TCGA-5T-A9QA-01A-11D-A41E-01 BMS1,RET,RASGEF1A,ZNF33B,CSGALNAC

15、T2chr10 106384615 106473355 TCGA-5T-A9QA-01A-11D-A41E-01 SORCS3chr10 106478366 107298256 TCGA-5T-A9QA-01A-11D-A41E-01 SORCS3chr10 117457285 117457859 TCGA-5T-A9QA-01A-11D-A41E-01 ATRNL1chr11 68990173 69277078 TCGA-5T-A9QA-01A-11D-A41E-01 MYEOVchr11 76378708 76926535 TCGA-5T-A9QA-01A-11D-A41E-01 LRRC32,B3GNT6,OMP,TSKU,MYO7