整理ChIP-seq / CUT & Tag 分析時(shí)用到的工具诵闭。本文只對(duì)使用的工具用法進(jìn)行簡(jiǎn)單介紹怎披。
一般而言誉券，我們獲得了MACS的peak calling結(jié)果時(shí)朦蕴，我們會(huì)想探究不同處理的peak calling結(jié)果是否會(huì)有重疊篮条。bedtools intersect就可以有效地解決這種問(wèn)題弟头。以下使用bedtools官方提供的例子進(jìn)行展示
一比一
一般情況，使用bedtools intersect尋找兩個(gè)peaks文件的交集可以如下操作
$ cat A.bed
chr1 10 20
chr1 30 40
$ cat B.bed
chr1 15 20
$ bedtools intersect -a A.bed -b B.bed
chr1 15 20
-a：查詢文件涉茧，可以是BAM/BED/GFF/VCF格式赴恨，A中的每一行都會(huì)與B進(jìn)行比較以找出重疊的部分
-b：參考文件，依照-b輸入的文件在A里面找重疊的部分降瞳，可以是BAM/BED/GFF/VCF格式嘱支。-b可以輸入多個(gè)文件，用空格分開(kāi)即可≌跫ⅲ現(xiàn)在可以用通配符（*）識(shí)別多個(gè)輸入文件。
一比多
$ cat query.bed
chr1 1 20
chr1 40 45
chr1 70 90
chr1 105 120
chr2 1 20
chr2 40 45
chr2 70 90
chr2 105 120
chr3 1 20
chr3 40 45
chr3 70 90
chr3 105 120
chr3 150 200
chr4 10 20
$ cat d1.bed
chr1 5 25
chr1 65 75
chr1 95 100
chr2 5 25
chr2 65 75
chr2 95 100
chr3 5 25
chr3 65 75
chr3 95 100
$ cat d2.bed
chr1 40 50
chr1 110 125
chr2 40 50
chr2 110 125
chr3 40 50
chr3 110 125
$ cat d3.bed
chr1 85 115
chr2 85 115
chr3 85 115
$ bedtools intersect -wa -wb \
-a query.bed \
-b d1.bed d2.bed d3.bed \
-sorted \
-filenames
chr1 1 20 d1.bed chr1 5 25
chr1 40 45 d2.bed chr1 40 50
chr1 70 90 d1.bed chr1 65 75
chr1 70 90 d3.bed chr1 85 115
chr1 105 120 d2.bed chr1 110 125
chr1 105 120 d3.bed chr1 85 115
chr2 1 20 d1.bed chr2 5 25
chr2 40 45 d2.bed chr2 40 50
chr2 70 90 d1.bed chr2 65 75
chr2 70 90 d3.bed chr2 85 115
chr2 105 120 d2.bed chr2 110 125
chr2 105 120 d3.bed chr2 85 115
chr3 1 20 d1.bed chr3 5 25
chr3 40 45 d2.bed chr3 40 50
chr3 70 90 d1.bed chr3 65 75
chr3 70 90 d3.bed chr3 85 115
chr3 105 120 d2.bed chr3 110 125
chr3 105 120 d3.bed chr3 85 115
-wa：在交集中輸出A的完整區(qū)間沛膳，而不僅是重疊的部分
-wb：在交集中輸出B的完整區(qū)間扔枫，而不僅是重疊的部分
-sorted：該參數(shù)表明輸入的文件是預(yù)先按照染色體和染色體坐標(biāo)排序過(guò)，可以提高bedtools運(yùn)行的效率锹安。
-filenames：標(biāo)明重疊區(qū)域是源于和哪個(gè)B文件發(fā)生重疊短荐。
輸出A unique的區(qū)域
$ bedtools intersect -wa -wb \
-a query.bed \
-b d1.bed d2.bed d3.bed \
-sorted \
-v
chr3 150 200
chr4 10 20
-v：輸出A unique的區(qū)域
定義重疊閾值
這里定義兩個(gè)區(qū)域100%一樣才算重疊
$ bedtools intersect -wa -wb \
-a query.bed \
-b d1.bed d2.bed d3.bed \
-sorted \
-names d1 d2 d3
-f 1.0
chr1 40 45 d2 chr1 40 50
chr2 40 45 d2 chr2 40 50
chr3 40 45 d2 chr3 40 50
-f：指定兩個(gè)區(qū)域最小重疊多少bp才算與A發(fā)生重疊，默認(rèn)1bp（1E9）
-F：指定兩個(gè)區(qū)域最小重疊多少bp才算與B發(fā)生重疊叹哭，默認(rèn)1bp（1E9）
輸出重疊區(qū)域的長(zhǎng)度
$ cat A.bed
chr1 10 20
chr1 30 40
$ cat B.bed
chr1 15 20
chr1 18 25
# -wa
$ bedtools intersect -a A.bed -b B.bed -wo
chr1 10 20 chr1 15 20 5
chr1 10 20 chr1 18 25 2
# -wao
$ 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 0
-wo：輸出A忍宋、B重疊區(qū)域重疊的bp
-wao：輸出A、B重疊區(qū)域重疊的bp风罩，假如B中不存在A的這個(gè)區(qū)域則輸出0糠排。
合并輸出結(jié)果中重疊的區(qū)域
有時(shí)候A與B的交集區(qū)域之間可能也有重疊，這一點(diǎn)可以通過(guò)-u參數(shù)解決超升。
$ cat A.bed
chr1 10 20
$ cat B.bed
chr1 15 20
chr1 17 22
# without `-u`
$ bedtools intersect -a A.bed -b B.bed
chr1 15 20
chr1 17 20
# with `-u`
$ bedtools intersect -a A.bed -b B.bed -u
chr1 10 20
輸出A與B overlap 的次數(shù)
如果我們只是想知道A文件中的區(qū)域與B文件的區(qū)域發(fā)生了多少次重疊入宦，可以使用-c參數(shù)，或者使用-C參數(shù)分別輸出不同來(lái)源的overlaps
cat A.bed
chr1 10 20
chr1 30 40
$ cat B.bed
chr1 15 20
chr1 18 25
$ cat C.bed
chr1 16 21
chr1 19 26
# -c 不管來(lái)源
$ bedtools intersect -a A.bed -b B.bed C.bed -c
chr1 10 20 4
chr1 30 40 0
# -C 考慮來(lái)源
$ bedtools intersect -a A.bed -b B.bed C.bed -C -filenames
chr1 10 20 B.bed 2
chr1 10 20 C.bed 2
chr1 30 40 B.bed 0
chr1 30 40 C.bed 0
考慮正負(fù)鏈
如果在找交集的過(guò)程中需要考慮鏈特異性的時(shí)候室琢，可以使用-s參數(shù)乾闰；
或者使用-S參數(shù)，強(qiáng)制不考慮正負(fù)鏈盈滴，而只考慮基因組區(qū)域是否重疊
$ cat A.bed
chr1 100 200 a1 100 +
$ cat B.bed
chr1 130 201 b1 100 -
chr1 132 203 b2 100 +
# 考慮正負(fù)鏈
$ bedtools intersect -a A.bed -b B.bed -wa -wb -s
chr1 100 200 a1 100 + chr1 132 203 b2 100 +
# 不考慮正負(fù)鏈涯肩，只考慮基因組區(qū)域
$ bedtools intersect -a A.bed -b B.bed -wa -wb -S
chr1 100 200 a1 100 + chr1 130 201 b1 100 -
對(duì)bedtools intersect的介紹到此結(jié)束，bedtools的其他命令有機(jī)會(huì)再做介紹巢钓。
ref
bedtools intersect docs: https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html
完病苗。
