确定了与表型显著的SNP位点后,通常会选择在显著性位点所在位置前后,各一定的距离内,确定候选区间,进行候选基因挖掘,而这个距离如何确定?
9.5 GWAS显著SNP筛选及曼哈顿图绘制 - 简书 (jianshu.com)
下文都以显著性位点:Chr1__272173254为例
方法一:主观设置
当测序密度较低时,基因组的覆盖度不够,得到的标记数据过少,标记之间的距离太大,无法构成LD block,这时可以分析师主观设定一个距离,如100k或更大,需要根据区间内基因的数目进行调整。当时这种方法的结果应该是最粗糙的。
此时候选区间就为:Chr1__272073254:Chr1__272273254
方法二:全基因组范围内的LD decay
通常我们会通过LD衰减距离来确定候选基因的筛选范围,即LD decay到r2等于0.6时所对应的距离,将GWAS超过阈值的SNP位点,即信号位点,前后各延伸这个距离作为候选区间,对这个区间内所有的基因进行注释,挑选候选基因。
10.GWAS:LD decay(LD衰减)—— PopLDdecay - 简书 (jianshu.com)
假设此时LD decay到r2等于0.6时距离为50k:
此时的候选区间为:Chr1__272123254:Chr1__272223254
这种方法是现在使用最广泛的,但是同时存在一个问题,基因组上的LD衰减速度并不是一致的,当时用全基因组的LD decay时,部分区域会因为连锁强度低于全基因组而漏掉一些关键基因,部分区域会因为连锁过强,导致候选基因增多,加大了后期的工作量。
方法三:部分区域内的LD decay
为解决方法二的缺陷,可以采取方法三的策略,即在方法二确定的候选区间内,计算区域内所有SNP的Pairwise correlations (r2),将r2大于0.6的block作为候选区间。
plink关于Pairwise correlations (r2)的官方说明书:
https://zzz.bwh.harvard.edu/plink/ld.shtml
$ plink --noweb --bfile root.id \
--chr 1 --from-bp 272123254 --to-bp 272223254 \
--r2 --out candidate_1.out \
--chr-set 27 \
--allow-extra-chr
--bfile:二进制bam文件的前缀
--out:输出文件的前缀
--chr-set:默认人类染色体,如果超出需要设置
--allow-extra-chr:允许有其他染色体,如scaffold等
# 查看结果文件
$ less candidate_1.out.ld
CHR_A BP_A SNP_A CHR_B BP_B SNP_B R2
1 272167626 Chr1__272167626 1 272173254 Chr1__272173254 0.357864
1 272167657 Chr1__272167657 1 272167710 Chr1__272167710 0.880447
1 272167657 Chr1__272167657 1 272173441 Chr1__272173441 0.525062
结果文件为指定区间内,所有SNP两两之间的R2。可以看到Chr1__272167657:Chr1__272167710的R2大宇0.6,因此候选区间即为Chr1__272167657:Chr1__272167710。
引用转载请注明出处,如有错误敬请指出。
网友评论