bedtools 用法大全

原文：https://cloud.tencent.com/developer/article/1078324

 

前言：

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