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

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

《bedtools用法大全一文就够吧.docx》由会员分享，可在线阅读，更多相关《bedtools用法大全一文就够吧.docx（7页珍藏版）》请在冰豆网上搜索。

bedtools用法大全一文就够吧

bedtools用法大全（一文就够吧）

bedtools等工具号称是可以代替普通的生物信息学数据处理工程师的！

我这里用一个专题来讲解它的用法，其实它能实现的需求，我们写脚本都是可以做的，而且我强烈建议正在学编程的小朋友模仿它的各种功能来增强自己的脚本功力。

BEDTools是可用于genomicfeatures的比较，相关操作及进行注释的工具。

而genomicfeatures通常使用BrowserExtensibleData（BED）或者GeneralFeatureFormat（GFF）文件表示，用UCSCGenomeBrowser进行可视化比较。

bedtools总共有二三十个工具/命令来处理基因组数据。

比较典型而且常用的功能举例如下：

格式转换，bam转bed（bamToBed），

bed转其他格式（bedToBam，bedToIgv）；

对基因组坐标的逻辑运算，包括：

交集（intersectBed，windowBed），”邻集“（closestBed），

补集（complementBed），并集（mergeBed），差集（subtractBed）;

计算覆盖度（coverage）（coverageBed，genomeCoverageBed）；

好，言归正传，bedtools是非常多的工具的合集，有瑞士军刀的美誉。

直接下载二进制版本软件就可以调用全路径来使用，或者把整个文件夹添加到环境变量。

首发于生信技能树：

genomecov我们先看第一个功能，把alignment的结果文件转为bedgraph格式文件。

不过这个功能用处不是很大。

参考：

http:

//bedtools.readthedocs.io/en/latest/content/tools/genomecov.html要实现这个功能，需要用bedtools的genomecov小命令，有两种方法可以调用：

bedtoolsgenomecov[OPTIONS][-i|-ibam]-g（iff.-i）

genomeCoverageBed[OPTIONS][-i|-ibam]-g（iff.-i）

这个命令本身并不是设计来做格式转换的，bam2bedgraph也只是其中的一个小功能而已，需要加上-bg参数，就可以ReportdepthinBedGraphformat.Fordetails,see:

http:

//genome.ucsc.edu/goldenPath/help/bedgraph.html大家观摩我下面给出的测试例子，就明白该功能如何使用了bedtoolsgenomecov-bg-iE001-H3K4me1.tagAlign-gmygenome.txt>E001-H3K4me1.bedGraph

bedtoolsgenomecov-bg-iE001-Input.tagAlign-gmygenome.txt>E001-Input.bedGraph

nohupbedtoolsgenomecov-bg-ibamBAF180_CT10.unique.sorted.bam>BAF180_CT10.bedGraph&

nohupbedtoolsgenomecov-bg-ibamBAF180_CT22AM.unique.sorted.bam>BAF180_CT22AM.bedGraph&

nohupbedtoolsgenomecov-bg-ibamBAF180_CT22.unique.sorted.bam>BAF180_CT22.bedGraph&

nohupbedtoolsgenomecov-bg-ibaminputCT10sonication.unique.sorted.bam>inputCT10sonication.bedGraph&

首先alignment的文件必须是sort的，然后如果是bed格式的比对文件，用-i参数来指定输入文件，需要加入参考基因组的染色体大小记录文件（mygenome.txt），如果是bam格式的比对文件，用-ibam指定输入文件，而且不需要参考基因组的染色体大小记录文件。

第二个功能multicov然后看看第二个功能,对RNA-seq的比对文件中的比对到各个基因的reads进行计数。

参考:

http:

//www.bio-info-这个小命令#例子：

bedtoolsmulticov-bamsaln1.bamaln2.bamaln3.bam-bedivls-of-interest.bed

#ivls-of-interest.bed这个文件是必须的，可能需要自己制作，其实用gtf文件也可以的，如下：

chr1010000ivl1

chr11000020000ivl2

chr12000030000ivl3

chr13000040000ivl4

输出结果前三列是坐标，第四列是基因名，跟我们的bed文件一样，只是最后三列是三个样本的计数，是添加上来的！

chr1010000ivl110022340

chr11000020000ivl212332451000

chr12000030000ivl321323322034

chr13000040000ivl433576540

可以看到，它实现的需求，跟htseq这个软件差不多。

各种区别，大家可以自己取探索。

第三个功能getfasta接着第三个功能，根据坐标区域来从基因组里面提取fasta序列参考：

#BED/GFF/VCF+reference-->fastabedtoolsgetfasta-fi~/biosoft/bowtie/hg19_index/hg19.fa-bed../macs14_results/highQuality_summits.bed-fohighQuality.fa

bedtoolsgetfasta-fi~/biosoft/bowtie/hg19_index/hg19.fa-bed../macs14_results/highQuality_peaks.bed-fohighQuality.fa

我的例子脚本里面用的是bed格式来记录坐标区域，参考基因组用-fi参数指定具体位置，输出的fasta序列文件用-fo参数指定第四个功能区域注释intersect首发于生信技能树论坛：

~/reference/gtf/gencode/protein_coding.hg19.position

chr16909170008OR4F5

chr1367640368634OR4F29

chr1621096622034OR4F16

chr1859308879961SAMD11

chr1879584894689NOC2L

chr1895967901095KLHL17

chr1901877911245PLEKHN1

chr1910584917473PERM1

chr1934342935552HES4

chr1936518949921ISG15

然后要把我们下载的CNV文本文件，转为bed格式的，就是把列的顺序调换一下：

headFeatures.bed

chr1321861095674710TCGA-3C-AAAU-10A-01D-A41E-01532250.0055

chr19567651195676518TCGA-3C-AAAU-10A-01D-A41E-012-1.6636

chr195680124167057183TCGA-3C-AAAU-10A-01D-A41E-01248860.0053

chr1167057495167059336TCGA-3C-AAAU-10A-01D-A41E-013-1.0999

chr1167059760181602002TCGA-3C-AAAU-10A-01D-A41E-019213-8e-04

chr1181603120181609567TCGA-3C-AAAU-10A-01D-A41E-016-1.2009

chr1181610685201473647TCGA-3C-AAAU-10A-01D-A41E-01120020.0055

chr1201474400201474544TCGA-3C-AAAU-10A-01D-A41E-012-1.4235

chr1201475220247813706TCGA-3C-AAAU-10A-01D-A41E-0129781-4e-04

命令很简单，如下：

bedtoolsintersect-aFeatures.bed-b~/reference/gtf/gencode/protein_coding.hg19.position\

-wa-wb|bedtoolsgroupby-i--g1-4-c10-ocollapse

注释结果如下：

可以看到，每个CNV片段都注释到了对应的基因，有些特别大的片段，会被注释到非常多的基因。

chr84258492442783715TCGA-5T-A9QA-01A-11D-A41E-01CHRNB3,CHRNA6,THAP1,RNF170,HOOK3

chr84278972842793594TCGA-5T-A9QA-01A-11D-A41E-01HOOK3

chr84279795742933372TCGA-5T-A9QA-01A-11D-A41E-01RP11-598P20.5,FNTA,HOOK3

chr87095267370964372TCGA-5T-A9QA-01A-11D-A41E-01PRDM14

chr104294797043833200TCGA-5T-A9QA-01A-11D-A41E-01BMS1,RET,RASGEF1A,ZNF33B,CSGALNACT2

chr10106384615106473355TCGA-5T-A9QA-01A-11D-A41E-01SORCS3

chr10106478366107298256TCGA-5T-A9QA-01A-11D-A41E-01SORCS3

chr10117457285117457859TCGA-5T-A9QA-01A-11D-A41E-01ATRNL1

chr116899017369277078TCGA-5T-A9QA-01A-11D-A41E-01MYEOV

chr117637870876926535TCGA-5T-A9QA-01A-11D-A41E-01LRRC32,B3GNT6,OMP,TSKU,MYO7