如何使用bedtools處理Rang數(shù)據(jù)
什么是Range數(shù)據(jù)
參考基因組表示的是一種坐標(biāo)系統(tǒng)苞也，比如說某一個(gè)物種基因組大小為100bp，那么他參考基因組就可以表示為[1,100], 之后就可以用任意[x,y]表示這條參考基因組上的位置粘秆，這就是一種范圍信息如迟，X-Y這段區(qū)域可能是外顯子，也可能是內(nèi)含子，可能是編碼區(qū)殷勘，也可能是基因間區(qū)此再，也有可能是一個(gè)測序結(jié)果。
因此Range數(shù)據(jù)是生信數(shù)據(jù)比較常見的存放形式玲销，比如說BED/BAM/BCF/和GFF/BFF/SAM/VCF/输拇，前者以0為始，后者以1為始贤斜。
為了操作這種Range數(shù)據(jù)策吠，Bioconductor在R語言中定義了兩個(gè)重要的對象，IRange和GenomicRanges瘩绒，后者僅存放'start','end','width'是后者的基礎(chǔ)猴抹。后者才能真正存放基因組Range數(shù)據(jù)。
這一篇不介紹如何在R語言操作Range數(shù)據(jù)草讶，而是介紹bedtools這款號稱基因組Range數(shù)據(jù)分析的瑞士軍當(dāng)洽糟，當(dāng)時(shí)的口號是一款取代10個(gè)生信分析師的工具。
Bedtools能夠?qū)蚪MRange數(shù)據(jù)進(jìn)行交堕战，并，補(bǔ)拍霜，計(jì)數(shù)等簡單操作嘱丢，也能和Unix命令行結(jié)合起來完成更加復(fù)雜的任務(wù)。
bed格式簡介
在正式介紹bedtools之前祠饺，需要先介紹一下BED格式越驻。根據(jù)USCSC基因組瀏覽器的描述，BED格式能夠非常簡潔的表示基因組特征和注釋道偷，盡管BED格式描述中定義了12列缀旁，但是僅僅只有3列必須，因此BED格式按照列數(shù)繼續(xù)細(xì)分為BED3,BED4,BED5,BED6,BED12朋鞍。
BED12定義的12列分別為：chrom, start, end, name(BED代表的特征名),score(范圍為0~1000亭饵，可以是pvalue, 或者是字符串,如"up"), strand(正負(fù)鏈), thickstart, thickednd(額外著色位置, 比如說表示外顯子), itemRgb(RGB顏色,如255,0,0), blockCount(區(qū)塊數(shù)量, 如外顯子), blockSizes(由逗號隔開的區(qū)塊大小), blockStarts(由逗號隔開的區(qū)塊起始位點(diǎn))锭魔。
知道了BED12后，就可以對BED的細(xì)分格式進(jìn)行舉例說明
- BED3：
chr1 11873 14409 - BED4:
chr1 11873 14409 uc001aaa.3 - BED5:
chr1 11873 14409 uc001aaa.3 0 - BED6:
chr1 11873 14409 uc001aaa.3 0 + - BED12:
chr1 11873 14409 uc001aaa.3 0 + 11873 12000 123,123,123 3 354,109,1189, 0,739,1347,
除了官方的BED定義外懊渡，BEDtools定義了BEDPE用來表示基因組不連續(xù)的特征，比如說結(jié)構(gòu)變異或者雙端測序的reads军拟。在定義中10列是必須的剃执，為chrom11, start1, end1, chrom2, start2, end2, name, score, strand1, strand2。 這之后可以增加任意多的其他列懈息。
其他BEDtools支持的格式說明：
- GFF3: https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md
- SAM/BAM: http://samtools.github.io/hts-specs/SAMv1.pdf
- VCF: http://samtools.github.io/hts-specs/VCFv4.2.pdf
Bedtools工具介紹
bedtools的功能非常強(qiáng)大肾档，試圖解決你所遇到的所有和基因組位置運(yùn)算的問題以及周邊問題：基因組運(yùn)算，多文件比較， PE數(shù)據(jù)操作怒见，格式轉(zhuǎn)換戒祠， Fasta數(shù)據(jù)操作， BAM工具速种， 統(tǒng)計(jì)學(xué)相關(guān)工具姜盈，其他小工具
其中最重要的選項(xiàng)是--help，一個(gè)強(qiáng)大的工具提供了許多參數(shù)配阵，需要勤讀幫助文檔馏颂。bedtools的官方文檔寫的非常優(yōu)秀，絕大部分工具都以圖解的方式形象的說明了每個(gè)參數(shù)的可能效果棋傍。因此我寫這篇文章的目的就是讓迫使自己去熟悉所有的工具而已救拉。
核心：基因組運(yùn)算
所謂的基因組運(yùn)算，就是看看看自己手頭拿到的區(qū)域和你感興趣的區(qū)域的關(guān)系如何瘫拣。bedtools提供了如下工具做一系列你想到或者你想不到的事情亿絮。
集合運(yùn)算：
- intersect: 交集，也就是看兩個(gè)區(qū)域的重疊
- window: 和intersect類似也是求交集麸拄，但是還會看上下游一定區(qū)間內(nèi)是否有重疊派昧。 比如說看某個(gè)基因的上下游是否有peak
- complement: 補(bǔ)集，根據(jù)一個(gè)已有區(qū)域得到另一個(gè)不重疊的區(qū)域
- substract: 差集拢切， 求一個(gè)區(qū)域減去另一個(gè)區(qū)域后的結(jié)果
- merge: 合集蒂萎，將相鄰的區(qū)間合并成一個(gè)
- cluster: 類似于merge，但是不做合并
區(qū)間統(tǒng)計(jì)：
- closest: 找到和目標(biāo)最近的特征區(qū)域淮椰。比如說ChIP-seq得到的peaks可以根據(jù)“距離其最近的基因就是他調(diào)控的基因”的假設(shè)進(jìn)行注釋五慈。
- coverage: 計(jì)算某個(gè)區(qū)間的覆蓋情況
- genomecov: 計(jì)算全部基因組的覆蓋情況
- map: 在某個(gè)區(qū)間的應(yīng)用其他函數(shù)， 比如說求和(sum), 均值(mean)
- annotate: 從其他一系列文件中統(tǒng)計(jì)給定區(qū)間的覆蓋情況
區(qū)間工具：
- flank: 根據(jù)已有特征區(qū)間主穗，得到兩翼位置的新區(qū)間
- slop: 根據(jù)已有特征區(qū)間泻拦，向外衍生
- shift: 特征區(qū)域整體偏移一定位置
- shuffle: 從基因組區(qū)間中隨機(jī)選擇某些位置
- sample: 從已有的區(qū)域隨機(jī)挑選
- makewindows: 從基因組上創(chuàng)建等長區(qū)間
其他
bedtools的核心工具就是上面幾個(gè)，剩下的都算是小輪子忽媒，解決了你手上輪子不夠多的煩惱争拐。
Fasta相關(guān)：
- getfasta: 根據(jù)提供的區(qū)間提取序列
- maskfasta: 根據(jù)提供的區(qū)間對序列進(jìn)行遮蓋
- nuc: 統(tǒng)計(jì)某個(gè)區(qū)間的堿基含量
BAM格式相關(guān)：
- 格式轉(zhuǎn)換：bamtobed, bedtobam, bamtofastq,bedpetobam,bed12tobed6
- multicov: 計(jì)算某些位點(diǎn)上不同BAM文件覆蓋
- tag: 計(jì)算一個(gè)BAM在不同位置上覆蓋
PE文件操作：
- pairtobed: 專門用于看其他格式在BEDPE格式上重疊情況
- pairtopair: 比較兩個(gè)BEDPE格式的重疊情況
統(tǒng)計(jì)學(xué)工具：主要是用以不同的統(tǒng)計(jì)學(xué)方法來衡量兩個(gè)區(qū)間的相似度，有三種： jaccard, reldist, fisher
除了以上猾浦，還有一些更加有趣的小工具陆错，比如說igv可以創(chuàng)建IGV自動(dòng)截圖的運(yùn)行腳本，links可以構(gòu)建能在UCSC基因組上打開的鏈接等金赦。
