响水做网站的,郯城县网站建设,网站需求分,餐饮加盟网站模板大家好#xff0c;我是邓飞。
今天星球的小伙伴问了一个问题#xff1a; 我现在在做GWAS分析#xff0c;现在已经找到性状关联的SNP位点#xff0c;下一步我如何根据position 找到基因呢#xff1f; 关于基因注释#xff0c;之前写过一些博客#xff0c;可以用到的软件…大家好我是邓飞。
今天星球的小伙伴问了一个问题 我现在在做GWAS分析现在已经找到性状关联的SNP位点下一步我如何根据position 找到基因呢 关于基因注释之前写过一些博客可以用到的软件有ANNOVAR、Bedtools今天回答了这个问题感觉excel也可以做基因注释了。
下面对我的回答进行进一步的阐述。
1. GWAS分析
GWAS分析之前写过一个Cookbook包括方方面面的内容了如果是小白推荐一遍看配套的视频一遍敲代码学习
录制了配套的视频教程前面的数据下载、软件安装、环境配置等相关视频免费观看后面的付费观看。对于想要快速学习的小白视频代码数据实操技术支持是比较快的一条路。 扫码查看视频教程
2显著SNP位点
做完GWAS分析后确定阈值然后小于阈值的位点都是显著性位点显著性位点最重要的两个信息 染色体 物理位置
有时候还包括snp的名称但是不是必填项只需要上面两个信息就可以知道显著snp在基因组上的位置了。
3配套基因组的gff文件
一般有基因组数据的物种有基因组的版本还有配套的gff或者gff3格式的文件文件的内容里面有 染色体 基因起始位置 基因终止位置 基因功能描述 ……
类似
4计算LD衰减距离
为何要计算LD衰减距离呢是为了知道显著snp代表的区间因为存在连锁所以衰减距离就是确定snp所代表的有效区间可以代表这个有效区间的变异。虽然snp不在基因上但是如果snp的衰减距离区间内比如上下50kb包含基因那也可以说明这个基因是显著影响性状的。
所以计算了LD衰减距离显著性snp的信息就变成了 染色体 有效区间起始位置 有效区间终止位置
5用excel注释显著性snp
我们把gff文件简化一下整理成excel格式 怎么用excel表格呢可以手动查看也可以编写一个函数。
话说上面的显著性位点一共就6个SNP手动搞就行了。
第一个snp区间是1染色体5-15这个区间有gene1 第二个snp区间是1染色体10-20这个区间有gene2不是完全包括但是有交集也算是 第三个snp没有基因 第四个snpgene4 第五个snp没有基因 第六个snp没有基因
所以这些snp一共注释的基因有gene1, gene2, gene4
6我有1000个显著性位点谢谢
如果位点很多这就需要用到软件了bedttols
「换到基因注释的领域看一下相关需求」 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
「基因区间文件」
共有5个基因区间文件分别是染色体开始位置终止位置基因名称。 cat gene_infor.ped chr1 1 14 gene1 chr1 17 19 gene2 chr1 45 82 gene3 chr1 88 93 gene4
2. 提取每个SNP上面的基因
「需求」 每个SNP一行 如果有基因在其区间放到右边如果没有基因返回空 如果一个SNP区间对应多个基因写成多行
代码 intersect交集 -a第一个位置信息表 -b第二个位置信息表 -loj以第一个为基准返回结果
结果可以看到第二个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
可以看到将没有匹配到基因的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
合并的结果 $ cat snp_infor_merge.ped chr1 5 20 chr1 30 40 chr1 80 90
然后和基因的信息进行合并
$ 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
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
结果可以看到SNP1有一个基因SNP2有2个基因SNP3没有基因……
6. 基因注释的不同玩法
把上面SNP的区间作为显著性SNP上下游的信息把基因的信息作为gff基因文件就可以进行基因注释了
上面的玩法都可以做。
「注意将gff格式整理为染色体开始位置结束位置基因信息
snp区间整理为染色体开始区间结束区间」
可以实现的功能 每个SNP区间内的基因 每个SNP全进内基因的个数 合并SNP区间内的基因 合并SNP区间内基因的个数
