显著SNP的基因注释教程！
大家好，我是邓飞。前几天推荐了一个快速学习GWAS的方法，还没看过的点击查看：所以GWAS学习看视频还是看代码？
GWAS分析中，我们用基因型数据（SNP）+表型数据，进行关联分析，得到显著性的SNP，这些SNP有染色体和物理位置，那么我们如何对SNP进行基因注释呢？即，我们如何得到显著SNP附近的基因。
一般一个物种，基因都已经注释过了，保存在gtf或者gff文件中，有物理位置，基因区间，基因的大体功能，我们可以用显著的SNP查找上下游附近的基因，这就是基因注释。
基因注释，有几步，比如确定显著SNP上下游多长，来查找基因，这就需要计算LD衰减距离：LD衰减图绘制--PopLDdecay，然后根据上下游去和gff文件合并，把区间内的基因找到，这就找到目标基因了。
下面，用一个例子，来介绍一下操作的方法：
下图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
「基因区间文件：」
共有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
2. 提取每个SNP上面的基因
「需求：」
- 每个SNP一行
- 如果有基因在其区间，放到右边，如果没有基因，返回空
- 如果一个SNP区间对应多个基因，写成多行
代码：
bedtools intersect -a snp_infor.ped -b gene_infor.ped -loj
- 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 .
结果可以看到，第二个SNP区间，对应两个基因，写成了两行。第三个SNP区间没有对应基因，用-1表示占位。共返回8行信息。
3. 返回有基因信息的SNP
如果不想要占位符，只想返回有基因的SNP信息，可以命令如下：
bedtools intersect -a snp_infor.ped -b gene_infor.ped -wa -wb
结果：
$ 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
可以看到，将没有匹配到基因的SNP删除了。
上面的信息中，有些SNP匹配到了多个基因，也就是基因是有重复的。
- 如果我们想看每个SNP匹配的基因情况，可以用上面的结果
- 如果我们想看一下共有多少无重复的基因匹配，就需要对SNP区间先合并
4. 合并SNP区间再匹配
合并命令：
bedtools merge -i snp_infor.ped >snp_infor_merge.ped
原始数据：
$ cat snp_infor.ped
chr1 5 15
chr1 10 20
chr1 30 40
chr1 80 90
chr1 110 120
chr1 115 125
合并的结果：
$ cat snp_infor_merge.ped
chr1 5 20
chr1 30 40
chr1 80 90
chr1 110 125
然后和基因的信息进行合并：
$ 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
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
结果可以看到，SNP1有一个基因，SNP2有2个基因，SNP3没有基因……
6. 基因注释的不同玩法
把上面SNP的区间，作为显著性SNP上下游的信息，把基因的信息作为gff基因文件，就可以进行基因注释了！
上面的玩法都可以做。
「注意，将gff格式整理为：染色体，开始位置，结束位置，基因信息；
snp区间整理为：染色体，开始区间，结束区间」
可以实现的功能：
- 每个SNP区间内的基因
- 每个SNP全进内基因的个数
- 合并SNP区间内的基因
- 合并SNP区间内基因的个数