一般操作中，我們使用最直接的方法，比如選擇上下游50k的作為區(qū)間。

這里介紹bedtools，他有個強大的功能，merge和intersect，把區(qū)間合并，然后和gff計算交集。下面介紹一下用法。

1. bedtools軟件介紹

軟件github地址：https://github.com/arq5x/bedtools2/

比較好的介紹文檔：https://www./post/cnposts/151/

總的來說，bedtools實用工具是一把瑞士軍刀，用于廣泛的基因組學(xué)分析任務(wù)。最廣泛使用的工具支持基因組算法：即基因組集合論。例如，bedtools允許人們以廣泛使用的基因組文件格式（如BAM、BED、GFF/GTF、VCF）從多個文件中交叉、合并、計數(shù)、補充和混洗基因組間隔。

雖然每個單獨的工具都設(shè)計用于執(zhí)行相對簡單的任務(wù)（例如，交叉兩個間隔文件），但通過在UNIX命令行上組合多個bedtools操作，可以進(jìn)行相當(dāng)復(fù)雜的分析。

2. bedtools軟件安裝

軟件下載最新版：https://github.com/arq5x/bedtools2/releases/tag/v2.30.0

 wget https://github.com/arq5x/bedtools2/releases/download/v2.30.0/bedtools-2.30.0.tar.gz

文件：bedtools-2.30.0.tar.gz

解壓，進(jìn)入文件夾，安裝：

tar zxvf bedtools-2.30.0.tar.gz
cd bedtools2
make

將bin文件夾的內(nèi)容，copy到~/bin一份。這樣可以直接調(diào)用了。

測試安裝成功：

$ bedtools
bedtools is a powerful toolset for genome arithmetic.

Version:   v2.30.0
About:     developed in the quinlanlab.org and by many contributors worldwide.
Docs:      http://bedtools.readthedocs.io/
Code:      https://github.com/arq5x/bedtools2
Mail:      https://groups.google.com/forum/#!forum/bedtools-discuss

Usage:     bedtools <subcommand> [options]

3. 交集常用方法

這里，我們構(gòu)建一個A.ped和B.ped。注意，分隔符為tab鍵。數(shù)據(jù)為三列：染色體，起始位置，終止位置。摸你的數(shù)據(jù)，和上圖的結(jié)構(gòu)一致。

「A.ped文件：」

$ cat A.ped
chr1 10 20
chr1 30 40
chr1 80 90

「B.ped文件：」

$ cat B.ped
chr1 1 14
chr1 17 19
chr1 45 82
chr1 88 93

3.1 默認(rèn)交集

命令：

bedtools intersect -a A.ped -b B.ped

• intersect，是交集
• -a，第一個ped文件
• -b，第二個ped文件

結(jié)果如下：

• 第一個重復(fù)區(qū)段是10~14
• 第二個重復(fù)區(qū)間是17~19
• 第三個重復(fù)區(qū)間是80~82
• 第四個重復(fù)區(qū)間是88~90

$ bedtools intersect -a A.ped -b B.ped
chr1 10 14
chr1 17 19
chr1 80 82
chr1 88 90

3.2 計算A中有B有交叉的區(qū)間

結(jié)果輸出時，加上參數(shù)-wa，返回的是A中的結(jié)果：

包含著染色體位置的兩個文件，分別記為A文件和B文件。求A文件中哪些染色體位置是與文件B中的染色體位置有overlap.

$ bedtools intersect -a A.ped -b B.ped -wa
chr1 10 20
chr1 10 20
chr1 80 90
chr1 80 90

3.3 計算B中有A有交叉的區(qū)間

包含著染色體位置的兩個文件，分別記為A文件和B文件。求A文件中染色體位置與文件B中染色體位置的交集，以及對應(yīng)的文件B中的染色體位置.

結(jié)果輸出時，加上參數(shù)-wb，返回B的結(jié)果：

$ bedtools intersect -a A.ped -b B.ped -wb
chr1 10 14 chr1 1 14
chr1 17 19 chr1 17 19
chr1 80 82 chr1 45 82
chr1 88 90 chr1 88 93

3.4 交集中，有的話返回B，沒有返回占位符

包含著染色體位置的兩個文件，分別記為A文件和B文件。求對于A文件的染色體位置是否與文件B中的染色體位置有交集。如果有交集，分別輸入A文件的染色體位置和B文件的染色體位置；如果沒有交集，輸入A文件的染色體位置并以'. -1 -1'補齊文件。

$ bedtools intersect -a A.ped -b B.ped -loj
chr1 10 20 chr1 1 14
chr1 10 20 chr1 17 19
chr1 30 40 . -1 -1
chr1 80 90 chr1 45 82
chr1 80 90 chr1 88 93

3.5 高階用法1

包含著染色體位置的兩個文件，分別記為A文件和B文件。對于A文件中染色體位置，如果和B文件中染色體位置有overlap,則輸出在A文件中染色體位置和在B文件中染色體位置，以及overlap的長度.

$ bedtools intersect -a A.ped -b B.ped -wo
chr1 10 20 chr1 1 14 4
chr1 10 20 chr1 17 19 2
chr1 80 90 chr1 45 82 2
chr1 80 90 chr1 88 93 2

3.6 高階用法2

包含著染色體位置的兩個文件，分別記為A文件和B文件。對于A文件中染色體位置，如果和B文件中染色體位置有overlap,則輸出在A文件中染色體位置和在B文件中染色體位置，以及overlap的長度；如果和B文件中染色體位置都沒有overlap,則用'. -1-1'補齊文件

$ bedtools intersect -a A.ped -b B.ped -wao
chr1 10 20 chr1 1 14 4
chr1 10 20 chr1 17 19 2
chr1 30 40 . -1 -1 0
chr1 80 90 chr1 45 82 2
chr1 80 90 chr1 88 93 2

3.7 高階用法3

包含著染色體位置的兩個文件，分別記為A文件和B文件。對于A文件中染色體位置，輸出在A文件中染色體位置和有多少B文件染色體位置與之有overlap.

$ bedtools intersect -a A.ped -b B.ped -c
chr1 10 20 2
chr1 30 40 0
chr1 80 90 2

3.8 高階用法4

包含著染色體位置的兩個文件，分別記為A文件和B文件。對于A文件中染色體位置，輸出在A文件中染色體位置和與B文件染色體位置至少有X%的overlap的記錄。

$ bedtools intersect -a A.ped -b B.ped -wa -wb -f 0.1
chr1 10 20 chr1 1 14
chr1 10 20 chr1 17 19
chr1 80 90 chr1 45 82
chr1 80 90 chr1 88 93

4. 顯著snp上下游確定

參考：https://www.jianshu.com/p/905bf1f3a798

第一種：每個顯著性位點，上下加100kb，為注釋區(qū)間。

直接當(dāng)測序密度較低時，基因組的覆蓋度不夠，得到的標(biāo)記數(shù)據(jù)過少，標(biāo)記之間的距離太大，無法構(gòu)成LD block，這時可以分析師主觀設(shè)定一個距離，如100k或更大，需要根據(jù)區(qū)間內(nèi)基因的數(shù)目進(jìn)行調(diào)整。當(dāng)時這種方法的結(jié)果應(yīng)該是最粗糙的。

第二種：每個顯著性位點，上下加LD衰減一般的距離為注釋區(qū)間 比如LD衰減圖是50kb，那就上下加50kb

第三種：根據(jù)顯著位點附件的LD decay

將大于0.6的block作為候選區(qū)間。

為了操作方便，我們使用LD衰減一般的距離為上下區(qū)間。

這里編寫一個腳本，自動添加上下游確定區(qū)間。生成三列的數(shù)據(jù)：染色體，開始區(qū)間，結(jié)束區(qū)間

$ head snp_re_qujian_sort.txt
1 44459430 44461430
1 44459431 44461431
1 114533676 114535676
1 114570742 114572742
1 114574836 114576836
1 114575784 114577784
1 114583582 114585582
1 115319771 115321771
3 2266311 2268311
3 2433000 2435000

5. 將區(qū)間去重合并

bedtools merge -i snp_re_qujian_sort.txt >snp_qujian.bed

$ head snp_qujian.bed
1 44459430 44461431
1 114533676 114535676
1 114570742 114572742
1 114574836 114577784
1 114583582 114585582
1 115319771 115321771
3 2266311 2268311
3 2433000 2435000
3 2457627 2459627
3 2460582 2462582

6. 將下載的gff3，進(jìn)行排序，然后與顯著性區(qū)間進(jìn)行合并

bedtools sort -chrThenSizeA -i aa.gff3.gz >t2.gff3

bedtools intersect -a t2.gff3 -b snp_qujian.bed -wa >re1.txt

這個re1.txt文件，即是注釋的區(qū)間，可以從中挑選基因信息和其它信息。