大家好，我是邓飞。
昨天推荐了一个GWAS培训，大家热情很高，没有看到的小伙伴快来了解一波：GWAS分析先做后学。
今天我把之前的GWAS教程更新了一般，工作量太大了，还没有搞完。我用Typora重新编辑了一下，界面美观多了，有增加了一些内容，明天下午做好之后会发个公众号，让大家领取，敬请期待呀！
今天介绍一下GWAS分析中注释的方法，我们知道，GWAS分析找到显著性SNP后，需要注释，才能找到候选的基因。
什么是注释呢？
我们知道，GWAS的依据是snp与控制性状的基因处于LD状态，所以，我们才能推断显著性的SNP附近的基因是影响性状的候选基因。那么这个附近如何确定呢？
1，最简单的方法，只看基因上的SNP显著位点，这个准确性是最高的，我们知道SNP有染色体和物理位置，而注释基因gff文件，也有基因的染色体和物理位置区间，我们就可以通过查看得到显著SNP所在的基因。当然，这种方法得到的候选基因最少，因为本来显著的SNP就比较少，而处在基因上的就更少了。这种方法也没有必要，因为显著SNP往往处于Block，强连锁，所以显著SNP附近的基因都有可能是候选基因。,
2，最直接的方法，把显著性的SNP，取LD的半衰期比如50kb，选择上下游50kb作为候选区间，然后和gff文件比对，处在这个区间内的基因，都是候选基因。这种方法，操作简单，但是不是每个染色体每个区段的LD都是50kb，因为我们算的是真个基因组的平均半衰期。
3，最准确的方法，根据每个显著SNP，计算他附近的LD值，选择一定的值把所有相关的SNP都找到，作为其上下游区间，比如下图把基因处于LDblock的snp都提取出来，和gff文件合并。但是这种方法过于复杂。
一般操作中，我们使用最直接的方法，比如选择上下游50k的作为区间。
这里介绍bedtools，他有个强大的功能，merge和intersect，把区间合并，然后和gff计算交集。下面介绍一下用法。
软件github地址：https://github.com/arq5x/bedtools2/
比较好的介绍文档：https://www.jieandze1314.com/post/cnposts/151/
总的来说，bedtools实用工具是一把瑞士军刀，用于广泛的基因组学分析任务。最广泛使用的工具支持基因组算法：即基因组集合论。例如，bedtools允许人们以广泛使用的基因组文件格式（如BAM、BED、GFF/GTF、VCF）从多个文件中交叉、合并、计数、补充和混洗基因组间隔。
虽然每个单独的工具都设计用于执行相对简单的任务（例如，交叉两个间隔文件），但通过在UNIX命令行上组合多个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-discuss
Usage: bedtools <subcommand> [options]
这里，我们构建一个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
命令：
bedtools intersect -a A.ped -b B.ped
结果如下：
$ bedtools intersect -a A.ped -b B.ped
chr1 10 14
chr1 17 19
chr1 80 82
chr1 88 90
结果输出时，加上参数-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
包含着染色体位置的两个文件，分别记为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
包含着染色体位置的两个文件，分别记为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
包含着染色体位置的两个文件，分别记为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
包含着染色体位置的两个文件，分别记为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
包含着染色体位置的两个文件，分别记为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
包含着染色体位置的两个文件，分别记为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
参考：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
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
这个re1.txt文件，即是注释的区间，可以从中挑选基因信息和其它信息。
联系客服
