计算与某一个位点连锁不平衡的D'值

  • 首选选择要计算的区域，这里我选择候选位点上下游1M的区域
grep -v "#" chr4.vcf |awk '{if($2>=68647209&&$2<=70647209){print$3}}'|grep -v "CM000096.5:69647209" >snp1.list
#提取目标区域的所有snp ID 并去掉我们用来与其他snp ID比较的位点
mkdir snp1  
echo "4 68647209 70647209"|tr " " "\t">snp1.position
bgzip chr4.vcf
tabix -p vcf chr4.vcf.gz
bcftools filter -R snp1.position chr4.vcf.gz>snp1/snp1.vcf    #这样会大大提高计算速度
  • 用plink计算D‘
cd snp1/
for i in `cat ../snp1.list`;do plink --vcf snp1.vcf --ld CM000096.5:69647209 $i;echo CM000096.5:69647209 $i `grep R-sq plink.log` >>snp1.ld;rm -f plink.log;done
©著作权归作者所有,转载或内容合作请联系作者
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

推荐阅读更多精彩内容