使用bedtools进行gwas基因注释

大家好，我是邓飞。

GWAS分析中，找到显著性SNP，需要注释，才能找到候选的基因。

什么是注释呢？
我们知道，GWAS的依据是snp与控制性状的基因处于LD状态，所以，我们才能推断显著性的SNP附近的基因是影响性状的候选基因。那么这个附近如何确定呢？

最简单的方法，只看基因上的SNP显著位点，这个准确性是最高的，我们知道SNP有染色体和物理位置，而注释基因gff文件，也有基因的染色体和物理位置区间，我们就可以通过查看得到显著SNP所在的基因。当然，这种方法得到的候选基因最少，因为本来显著的SNP就比较少，而处在基因上的就更少了。这种方法也没有必要，因为显著SNP往往处于Block，强连锁，所以显著SNP附近的基因都有可能是候选基因。
最直接的方法，把显著性的SNP，取LD的半衰期比如50kb，选择上下游50kb作为候选区间，然后和gff文件比对，处在这个区间内的基因，都是候选基因。这种方法，操作简单，但是不是每个染色体每个区段的LD都是50kb，因为我们算的是真个基因组的平均半衰期。
最准确的方法，根据每个显著SNP，计算他附近的LD值，选择一定的值把所有相关的SNP都找到，作为其上下游区间，比如下图把基因处于LDblock的snp都提取出来，和gff文件合并。但是这种方法过于复杂。

一般操作中，我们使用最直接的方法，比如选择上下游50k的作为区间。

这里介绍bedtools，他有个强大的功能，merge和intersect，把区间合并，然后和gff计算交集。下面介绍一下用法：

1. bedtools软件介绍

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

比较好的介绍文档：https://www.jieandze1314.com/post/cnposts/151/

总的来说，bedtools实用工具是一把瑞士军刀，用于广泛的基因组学分析任务。最广泛使用的工具支持基因组算法：即基因组集合论。例如，bedtools允许人们以广泛使用的基因组文件格式（如BAM、BED、GFF/GTF、VCF）从多个文件中交叉、合并、计数、补充和混洗基因组间隔。

虽然每个单独的工具都设计用于执行相对简单的任务（例如，交叉两个间隔文件），但通过在UNIX命令行上组合多个bedtools操作，可以进行相当复杂的分析。

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

解压，进入文件夹，安装：

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

将bin文件夹的内容，copy到~/bin一份。这样可以直接调用了。

测试安装成功：

$ 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-discussUsage:     bedtools <subcommand> [options]

3. 交集常用方法

这里，我们构建一个A.ped和B.ped。注意，分隔符为tab键。数据为三列：染色体，起始位置，终止位置。摸你的数据，和上图的结构一致。

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 默认交集

命令：

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

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

结果如下：

第一个重复区段是10~14
第二个重复区间是17~19
第三个重复区间是80~82
第四个重复区间是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有交叉的区间

结果输出时，加上参数-wa，返回的是A中的结果：

包含着染色体位置的两个文件，分别记为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有交叉的区间

包含着染色体位置的两个文件，分别记为A文件和B文件。求A文件中染色体位置与文件B中染色体位置的交集，以及对应的文件B中的染色体位置.

结果输出时，加上参数-wb，返回B的结果：

$ 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，为注释区间。

直接当测序密度较低时，基因组的覆盖度不够，得到的标记数据过少，标记之间的距离太大，无法构成LD block，这时可以分析师主观设定一个距离，如100k或更大，需要根据区间内基因的数目进行调整。当时这种方法的结果应该是最粗糙的。

第二种：每个显著性位点，上下加LD衰减一般的距离为注释区间
比如LD衰减图是50kb，那就上下加50kb

第三种：根据显著位点附件的LD decay

将大于0.6的block作为候选区间。

为了操作方便，我们使用LD衰减一般的距离为上下区间。

这里编写一个脚本，自动添加上下游确定区间。生成三列的数据：染色体，开始区间，结束区间

$ 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. 将区间去重合并

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，进行排序，然后与显著性区间进行合并

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

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

这个re1.txt文件，即是注释的区间，可以从中挑选基因信息和其它信息。