使用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

1
2

文件：
bedtools-2.30.0.tar.gz

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

tar zxvf bedtools-2.30.0.tar.gz
cd bedtools2
make
1
2
3

将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-discuss

Usage:     bedtools <subcommand> [options]

1
2
3
4
5
6
7
8
9
10
11

3. 交集常用方法

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

A.ped文件：

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

1
2
3
4
5

B.ped文件：

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

1
2
3
4
5
6

3.1 默认交集

命令：

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

1
2

• 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

1
2
3
4
5
6

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

1
2
3
4
5
6

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
1
2
3
4
5

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

1
2
3
4
5
6
7

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

1
2
3
4
5
6

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

1
2
3
4
5
6
7

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

1
2
3
4
5

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


1
2
3
4
5
6
7

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


1
2
3
4
5
6
7
8
9
10
11
12
13

5. 将区间去重合并

bedtools merge -i snp_re_qujian_sort.txt >snp_qujian.bed
1

$ 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

1
2
3
4
5
6
7
8
9
10
11
12

6. 将下载的gff3，进行排序，然后与显著性区间进行合并

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

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

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