大家好，我是鄧飛。
昨天推薦了一個(gè)GWAS培訓(xùn)，大家熱情很高，沒(méi)有看到的小伙伴快來(lái)了解一波：GWAS分析先做后學(xué)。
今天我把之前的GWAS教程更新了一般，工作量太大了，還沒(méi)有搞完。我用Typora重新編輯了一下，界面美觀多了，有增加了一些內(nèi)容，明天下午做好之后會(huì)發(fā)個(gè)公眾號(hào)，讓大家領(lǐng)取，敬請(qǐng)期待呀！
今天介紹一下GWAS分析中注釋的方法，我們知道，GWAS分析找到顯著性SNP后，需要注釋?zhuān)拍苷业胶蜻x的基因。
什么是注釋呢？
我們知道，GWAS的依據(jù)是snp與控制性狀的基因處于LD狀態(tài)，所以，我們才能推斷顯著性的SNP附近的基因是影響性狀的候選基因。那么這個(gè)附近如何確定呢？
1，最簡(jiǎn)單的方法，只看基因上的SNP顯著位點(diǎn)，這個(gè)準(zhǔn)確性是最高的，我們知道SNP有染色體和物理位置，而注釋基因gff文件，也有基因的染色體和物理位置區(qū)間，我們就可以通過(guò)查看得到顯著SNP所在的基因。當(dāng)然，這種方法得到的候選基因最少，因?yàn)楸緛?lái)顯著的SNP就比較少，而處在基因上的就更少了。這種方法也沒(méi)有必要，因?yàn)轱@著SNP往往處于Block，強(qiáng)連鎖，所以顯著SNP附近的基因都有可能是候選基因。,
2，最直接的方法，把顯著性的SNP，取LD的半衰期比如50kb，選擇上下游50kb作為候選區(qū)間，然后和gff文件比對(duì)，處在這個(gè)區(qū)間內(nèi)的基因，都是候選基因。這種方法，操作簡(jiǎn)單，但是不是每個(gè)染色體每個(gè)區(qū)段的LD都是50kb，因?yàn)槲覀兯愕氖钦鎮(zhèn)€基因組的平均半衰期。
3，最準(zhǔn)確的方法，根據(jù)每個(gè)顯著SNP，計(jì)算他附近的LD值，選擇一定的值把所有相關(guān)的SNP都找到，作為其上下游區(qū)間，比如下圖把基因處于LDblock的snp都提取出來(lái)，和gff文件合并。但是這種方法過(guò)于復(fù)雜。
一般操作中，我們使用最直接的方法，比如選擇上下游50k的作為區(qū)間。
這里介紹bedtools，他有個(gè)強(qiáng)大的功能，merge和intersect，把區(qū)間合并，然后和gff計(jì)算交集。下面介紹一下用法。
軟件github地址：https://github.com/arq5x/bedtools2/
比較好的介紹文檔：https://www.jieandze1314.com/post/cnposts/151/
總的來(lái)說(shuō)，bedtools實(shí)用工具是一把瑞士軍刀，用于廣泛的基因組學(xué)分析任務(wù)。最廣泛使用的工具支持基因組算法：即基因組集合論。例如，bedtools允許人們以廣泛使用的基因組文件格式（如BAM、BED、GFF/GTF、VCF）從多個(gè)文件中交叉、合并、計(jì)數(shù)、補(bǔ)充和混洗基因組間隔。
雖然每個(gè)單獨(dú)的工具都設(shè)計(jì)用于執(zhí)行相對(duì)簡(jiǎn)單的任務(wù)（例如，交叉兩個(gè)間隔文件），但通過(guò)在UNIX命令行上組合多個(gè)bedtools操作，可以進(jìn)行相當(dāng)復(fù)雜的分析。
軟件下載最新版：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)用了。
測(cè)試安裝成功：
$ 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]
這里，我們構(gòu)建一個(gè)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
命令：
bedtools intersect -a A.ped -b B.ped
結(jié)果如下：
$ bedtools intersect -a A.ped -b B.ped
chr1 10 14
chr1 17 19
chr1 80 82
chr1 88 90
結(jié)果輸出時(shí)，加上參數(shù)-wa，返回的是A中的結(jié)果：
包含著染色體位置的兩個(gè)文件，分別記為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
包含著染色體位置的兩個(gè)文件，分別記為A文件和B文件。求A文件中染色體位置與文件B中染色體位置的交集，以及對(duì)應(yīng)的文件B中的染色體位置.
結(jié)果輸出時(shí)，加上參數(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
包含著染色體位置的兩個(gè)文件，分別記為A文件和B文件。求對(duì)于A文件的染色體位置是否與文件B中的染色體位置有交集。如果有交集，分別輸入A文件的染色體位置和B文件的染色體位置；如果沒(méi)有交集，輸入A文件的染色體位置并以'. -1 -1'補(bǔ)齊文件。
$ 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
包含著染色體位置的兩個(gè)文件，分別記為A文件和B文件。對(duì)于A文件中染色體位置，如果和B文件中染色體位置有overlap,則輸出在A文件中染色體位置和在B文件中染色體位置，以及overlap的長(zhǎng)度.
$ 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
包含著染色體位置的兩個(gè)文件，分別記為A文件和B文件。對(duì)于A文件中染色體位置，如果和B文件中染色體位置有overlap,則輸出在A文件中染色體位置和在B文件中染色體位置，以及overlap的長(zhǎng)度；如果和B文件中染色體位置都沒(méi)有overlap,則用'. -1-1'補(bǔ)齊文件
$ 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
包含著染色體位置的兩個(gè)文件，分別記為A文件和B文件。對(duì)于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
包含著染色體位置的兩個(gè)文件，分別記為A文件和B文件。對(duì)于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
參考：https://www.jianshu.com/p/905bf1f3a798
第一種：每個(gè)顯著性位點(diǎn)，上下加100kb，為注釋區(qū)間。
直接當(dāng)測(cè)序密度較低時(shí)，基因組的覆蓋度不夠，得到的標(biāo)記數(shù)據(jù)過(guò)少，標(biāo)記之間的距離太大，無(wú)法構(gòu)成LD block，這時(shí)可以分析師主觀設(shè)定一個(gè)距離，如100k或更大，需要根據(jù)區(qū)間內(nèi)基因的數(shù)目進(jìn)行調(diào)整。當(dāng)時(shí)這種方法的結(jié)果應(yīng)該是最粗糙的。
第二種：每個(gè)顯著性位點(diǎn)，上下加LD衰減一般的距離為注釋區(qū)間 比如LD衰減圖是50kb，那就上下加50kb
第三種：根據(jù)顯著位點(diǎn)附件的LD decay
將大于0.6的block作為候選區(qū)間。
為了操作方便，我們使用LD衰減一般的距離為上下區(qū)間。
這里編寫(xiě)一個(gè)腳本，自動(dòng)添加上下游確定區(qū)間。生成三列的數(shù)據(jù)：染色體，開(kāi)始區(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
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
bedtools sort -chrThenSizeA -i aa.gff3.gz >t2.gff3
bedtools intersect -a t2.gff3 -b snp_qujian.bed -wa >re1.txt
這個(gè)re1.txt文件，即是注釋的區(qū)間，可以從中挑選基因信息和其它信息。
聯(lián)系客服
