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

(0)

相关推荐

  • Harvard FAS Informatics出品的ATAC

    Harvard FAS Informatics出品的ATAC-seq测序指南 github链接:harvardinformatics/ATAC-seq 参考文献:ATAC-seq: A Method ...

  • 【软件介绍】IGV软件的安装和基本介绍

    [软件介绍]IGV软件的安装和基本介绍 - 目录 1. IGV 下载与安装 2. Java 安装及环境变量设置 3. IGV 基本介绍 3.1 IGV 界面布局 3.2 IGV 结果界面 3.3 序列 ...

  • R语言版本的bedtools--bedtoolsr

    bedtools是一款非常强大的用于处理bed,vcf,gff等格式数据的工具,该软件由犹他大学的Quinlan实验室开发.但是目前bedtools主要提供的是在linux,unxi等操作系统环境下的 ...

  • pybedtools:对bedtools的封装和扩展

    bedtools是区间操作最常用的软件,pybedtools对其进行了封装,可以在python编程环境中灵活使用bedtools,而且进一步拓展出了很多有用的功能. 在pybedtools中, 核心是 ...

  • 英语数词用法大全(建议收藏)

      一.表示日期   1949年1949读作nineteen forty-nine 6月23日June 23rd 读作June (the) twenty-third或the twenty-third ...

  • 超全!各类设备原理,海量动图一文就够!

    这些展示设备原理的动图,一定有你没见过的. ▲ 火力发电流程原理 ▲ 核能发电流程原理 ▲ 水力发电流程原理 ▲ 光热发电原理 ▲ 垃圾发电原理 ▲ 蒸汽吸收式制冷原理 ▲ 尿素热解脱硝流程原理 ▲ ...

  • 函数 | COUNTIF用法大全

    COUNTIF作为最常用的函数之一,出现的也是最早的一批,在COUNTIFS没出现之前,他是计数的霸主,经过这么多年,被各位前辈开发出了很多用法,今天我们就来给大家盘上一盘. . 语法:COUNTIF ...

  • COUNTIF用法大全

    COUNTIF作为最常用的函数之一,出现的也是最早的一批,在COUNTIFS没出现之前,他是计数的霸主,经过这么多年,被各位前辈开发出了很多用法,今天我们就来给大家盘上一盘. . 语法:COUNTIF ...

  • 人工干预如何提高模型性能?看这文就够了!

    有一些行业对误报非常敏感,如金融行业,在对信用卡欺诈检测时,如果检测系统将用户的行为错误地分类为欺诈,这将对该金融机构的声誉产生负面影响[1].又如在医学领域对癌症诊断时,对假阳性反应是很敏感的[3] ...

  • 高中历史,如何用一句话搞定关键史实和核心概念?这一文就够了

    学习历史和语文一样,知识素养的积累需要一个长期的过程. 其实历史知识有其独特的规律,如果把相关知识归纳汇总,记忆起来也是很有效果的.今天我整理了中国史和世界史的关键史实和核心概念,简短精炼,方便记忆. ...

  • 给你总结全了:78个麻黄汤!(麻黄用法大全)

    麻黄汤是中医方剂的名称,是解表剂,麻黄汤具有发汗解表,宣肺平喘的功效,主治外感风寒表实证.恶寒发热.头身疼痛.无汗而喘.舌苔薄白.脉浮紧.临床上常用于治疗感冒.流行性感冒.急性支气管炎.支气管哮喘等属 ...

  • 【实用英文】你会说“sorry”吗?“sorry”花式用法大全!Different Ways to say “sorry”

    "sorry"是非常常用的一个词,几乎每个人每天都在用.它可以表示"抱歉:惭愧:惋惜:遗憾:难过"等种种意思. 大家最常说出口的是单字句: Sorry. 对不起 ...

  • 茶知识入门大全,这篇够份量!

    请识别下方二维码加我好友 一.中国 · 六大茶类 中国的六大茶类分别是:绿茶.红茶.白茶.青茶(乌龙茶).黄茶和黑茶. 整体而言,中国还是绿茶的天下,种类最繁多,毕竟江南自古善于自我营销,众才子佳人纷 ...