基因注释：区间范围匹配

大家好，我是邓飞。

今天，有老师问我一个问题，如果从一个区间匹配到另一个的区间范围，并找出来。我觉得比较有代表性，就写篇博客总结一下。

老师的需求如下：

图1是SNP的上下游区间，图2是基因的上下游区间，想以图1为标准，将区间内有基因的行放到右边。

换到基因注释的领域，看一下相关需求：

• 1，显著性的SNP位点，取上下游50k的位点，作为候选的区间
• 2，将候选区间有基因的，匹配到SNP的右边

处理注意：

• 1，显著SNP在上下游区间时，可能会有交叉，所以要先合并（merge）
• 2，匹配基因时，一个SNP区间可能会有多个基因

1. 数据描述

SNP区间文件：

这里，提取显著SNP的区间，提取三列信息：染色体，开始位置，结束位置：

共有6个SNP区间，其中第一个和第二个有重合，第五个和第六个有重合。

$ cat snp_infor.ped
chr1	5	15
chr1	10	20
chr1	30	40
chr1	80	90
chr1	110	120
chr1	115	125

1
2
3
4
5
6
7
8

基因区间文件：

共有5个基因区间文件，分别是：染色体，开始位置，终止位置，基因名称。

$ cat gene_infor.ped
chr1	1	14	gene1
chr1	17	19	gene2
chr1	45	82	gene3
chr1	88	93	gene4
chr1	100	105	gene5

1
2
3
4
5
6
7

2. 提取每个SNP上面的基因

需求：

• 每个SNP一行
• 如果有基因在其区间，放到右边，如果没有基因，返回空
• 如果一个SNP区间对应多个基因，写成多行

代码：

bedtools intersect -a snp_infor.ped -b gene_infor.ped  -loj

1
2

• intersect，交集
• -a，第一个位置信息表
• -b，第二个位置信息表
• -loj，以第一个为基准，返回结果

结果：

$ bedtools intersect -a snp_infor.ped -b gene_infor.ped  -loj
chr1	5	15	chr1	1	14	gene1
chr1	10	20	chr1	1	14	gene1
chr1	10	20	chr1	17	19	gene2
chr1	30	40	.	-1	-1	.
chr1	80	90	chr1	45	82	gene3
chr1	80	90	chr1	88	93	gene4
chr1	110	120	.	-1	-1	.
chr1	115	125	.	-1	-1	.

1
2
3
4
5
6
7
8
9
10

结果可以看到，第二个SNP区间，对应两个基因，写成了两行。第三个SNP区间没有对应基因，用-1表示占位。共返回8行信息。

3. 返回有基因信息的SNP

如果不想要占位符，只想返回有基因的SNP信息，可以命令如下：

bedtools intersect -a snp_infor.ped -b gene_infor.ped  -wa -wb

1
2

结果：

$ bedtools intersect -a snp_infor.ped -b gene_infor.ped  -wa -wb
chr1	5	15	chr1	1	14	gene1
chr1	10	20	chr1	1	14	gene1
chr1	10	20	chr1	17	19	gene2
chr1	80	90	chr1	45	82	gene3
chr1	80	90	chr1	88	93	gene4

1
2
3
4
5
6
7

可以看到，将没有匹配到基因的SNP删除了。

上面的信息中，有些SNP匹配到了多个基因，也就是基因是有重复的。

• 如果我们想看每个SNP匹配的基因情况，可以用上面的结果
• 如果我们想看一下共有多少无重复的基因匹配，就需要对SNP区间先合并

4. 合并SNP区间再匹配

合并命令：

bedtools merge -i snp_infor.ped >snp_infor_merge.ped

1
2

原始数据：

$ cat snp_infor.ped
chr1	5	15
chr1	10	20
chr1	30	40
chr1	80	90
chr1	110	120
chr1	115	125

1
2
3
4
5
6
7
8

合并的结果：

$ cat snp_infor_merge.ped
chr1	5	20
chr1	30	40
chr1	80	90
chr1	110	125

1
2
3
4
5
6

然后和基因的信息进行合并：

$ bedtools intersect -a snp_infor_merge.ped -b gene_infor.ped -wa -wb
chr1	5	20	chr1	1	14	gene1
chr1	5	20	chr1	17	19	gene2
chr1	80	90	chr1	45	82	gene3
chr1	80	90	chr1	88	93	gene4

1
2
3
4
5
6

5. 查看每个SNP区间基因的个数

结果可以用2中，统计一下个数，也可以用bedtools的-c参数：

$ bedtools intersect -a snp_infor.ped -b gene_infor.ped -c
chr1	5	15	1
chr1	10	20	2
chr1	30	40	0
chr1	80	90	2
chr1	110	120	0
chr1	115	125	0

1
2
3
4
5
6
7
8

结果可以看到，SNP1有一个基因，SNP2有2个基因，SNP3没有基因……

6. 基因注释的不同玩法

把上面SNP的区间，作为显著性SNP上下游的信息，把基因的信息作为gff基因文件，就可以进行基因注释了！

上面的玩法都可以做。

注意，将gff格式整理为：染色体，开始位置，结束位置，基因信息；snp区间整理为：染色体，开始区间，结束区间

可以实现的功能：

• 每个SNP区间内的基因
• 每个SNP全进内基因的个数
• 合并SNP区间内的基因
• 合并SNP区间内基因的个数