单倍型分析
理论上来说，单倍型的类别是2N，但由于LD的存在，会使得单倍型的类别远远的小于2N，但是如果当N的数目很大，即便是类别少，在整个群体中进行分类也会有很多种单倍型。这样的情况，可以对单倍型分成单倍型域，单体型域内几乎不发生重组

image.png

Hapview计算block和tagSNP

Hapview文件格式的要求

test.ped文件要求

image.png

该文件的第一列是家系的ID，如果是无关个体之间的研究，该列应为不重复的ID号；

第二列为个体的ID，做无关个体的研究，每个个体的编号不能重复；

第三列是父亲的ID，如果是无关个体的研究，该列为0；

第四列是母亲的ID，如果是无关个体的研究，该列为0；

第五列是个体的性别信息，其中1代表男性，2代表女性，如果研究个体的性别未知，该列用0表示；

第六列是个体的患病情况，0代表患病情况未知，1代表未患病， 2代表患病；

第七列之后，每两列代表个体一个位点的SNP信息，该软件不能识别A/T/C/G/N碱基，因此后续的SNP信息只能用数字表示， 其中1=A; 2=C; 3=G; 4=T; 缺失碱基用0表示；

sample.info文件

image.png

文件格式转换

##对vcf文件转换

###读入数据
rm(list=ls())
# choose.files()
library(vcfR)
setwd("D:\\mywork\\mywork\\16.小麦穗型世界和中国范围内变化GWAS\\1.321鲁非GWAS分析结果\\8.株型结果\\12.株高1A结果\\")
vcf_data <- read.vcfR("SNP_246_snp_plant_height_1A.chr1A.45791400-49094709.recode.vcf", verbose = FALSE )

##VCF文件包含GT:AD:GL，GT 代表的是genotype；AD代表的是allel depth；DP(Depth)为sample中该位点的覆盖度，是所支持的两个AD值（逗号前和逗号后）的加
#GL：三种基因型（RR RA AA）出现的可能性，R表示参考碱基，A表示变异碱基

vcf_data_gt <- extract.gt(vcf_data ,element="GT")


dim(vcf_data_gt)
vcf_data_fix <- vcf_data@fix
##对数据进行替换
##首先将缺失替换成0.5

vcf_data_gt[is.na(vcf_data_gt)] <- "0 0"
for (i in 1:9875) {
  vcf_data_gt[i,][vcf_data_gt[i,] == '0/0']<- paste(vcf_data_fix[i,4],vcf_data_fix[i,4],sep=" ")
  vcf_data_gt[i,][vcf_data_gt[i,] == '1/1']<- paste(vcf_data_fix[i,5],vcf_data_fix[i,5],sep=" ")
  ##将杂合也替换成0.5
  vcf_data_gt[i,][vcf_data_gt[i,] == '0/1']<- paste(vcf_data_fix[i,4],vcf_data_fix[i,5],sep=" ")
  vcf_data_gt[i,][vcf_data_gt[i,] == '1/0']<- paste(vcf_data_fix[i,5],vcf_data_fix[i,4],sep=" ")
  #
  
}

for (i in 1:9875) {
  vcf_data_gt[i,][vcf_data_gt[i,] == 'A A']<- "1 1"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'C C']<- "2 2"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'G G']<- "3 3"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'T T']<- "4 4"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'A C']<- "1 2"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'A G']<- "1 3"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'A T']<- "1 4"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'C A']<- "2 1"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'C G']<- "2 3"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'C T']<- "2 4"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'G A']<- "3 1"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'G C']<- "3 2"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'G T']<- "3 4"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'T A']<- "4 1"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'T C']<- "4 2"
  vcf_data_gt[i,][vcf_data_gt[i,] == 'T G']<- "4 3"
  
  
  
}

vcf_data_use<-t(vcf_data_gt)


write.csv(vcf_data_use,"vcf_data_use_HAPVIEW.csv")

直接将输出的数据按照格式要求整理好后，复制到他给的例子文件种，昨天进行这一步用自己将txt文件变成.ped文件一直出问题
Hapview如何输入数据如何进行参数的调整，可以参考下边连接
单倍型分析软件Haploview的导入格式及使用 - 云+社区 - 腾讯云 (tencent.com)
得到block绘图

image.png

存在的问题
（1）Hapview中计算block四种算法的区别还没研究透
（2）Hapview 存在java运行慢，SNP数目要求小，如果输入的数据多的话，导致软件直接卡死

用plink得到tagSNP

###首先得到目标区段的snp


sh get_gene_snp.sh steam_3D.txt chr3D.vcf

##将vcf文件修改成bed格式
plink --vcf steam4_peak_snp.chr2A.613558313-619417320.recode.vcf --recode --out test--const-fid --allow-extra-chr
plink --allow-extra-chr --noweb -file test--const-fid --out test1

##计算block
 plink --bfile test1 --blocks no-pheno-req --allow-extra-chr
##如果是bed格式就用--bfile
##如果是ped格式就用--file


##计算tagSNP
##首先准备一个block SNP文件信息，每一行是一个snp

##查看steam4_tagSNP
plink --bfile test1 --list-all --show-tags steam4_tagSNP.txt  --allow-extra-chr

[理学]动物数量性状单倍体型分析及其应用.ppt-原创力文档 (book118.com)
干货：如何利用1000 Genomes数据挑选tagSNP？ | 基因有限公司 (genecompany.com)
单倍型分析软件Haploview的导入格式及使用 - 云+社区 - 腾讯云 (tencent.com)
SNP连锁不平衡度量和单倍型构建 | Yang's Garden (yangli.name)
人类基因组的Phasing原理是什么？ - 知乎 (zhihu.com)
单倍体型分析方法及生物医学文献猪繁殖性状候选基因挖掘 - 豆丁网 (docin.com)
PLINK 连锁不平衡分析LD - 知乎 (zhihu.com)
(13条消息) haploview进行连锁不平衡分析_庐州月光的博客-CSDN博客

(13条消息) 使用shapeit进行单倍型分析_庐州月光的博客-CSDN博客

(13条消息) 使用Eagle2进行单倍型分析_庐州月光的博客-CSDN博客
基于RAINBOW的单倍型全基因组关联分析（haplotype-based GWAS）教程 - 知乎 (zhihu.com)

【翻译】RAINBOW：采用新型SNP-set方法的基于单倍型的全基因组关联分析【第二部分：材料和方法】 - Minerw - 博客园 (cnblogs.com)

用Beagle做基因型填充（Imputation） - 简书 (jianshu.com)
SR4R数据库：水稻4个SNP集的筛选及其应用 - 小xuo生 - 博客园 (cnblogs.com)

OmicShare Forum 广州基迪奥生物信息论坛 - Powered by Discuz!

基因组数据的定向和填充（Phasing and Imputation） - 简书 (jianshu.com)