目标基因单倍型分析——haploview软件使用及其数据整理

一，候选基因区间获得及对应区间VCF文件变异位点提取

1.tbtools提取x18参考基因组1号染色体基因组文件

2.转换为cds文件并进行blast，将berg参考基因组中得到的基因cds序列使用前面得到的1号染色体cds集合进行比对获得X18基因组中候选基因，使用X18注释文件获取候选基因的起始终止位置。

3.bcftoolst提取候选基因区间所有变异位点

# 激活环境，安装有bcftools,vcftools,plink,gemma分析软件
conda activate gwas

 bcftools filter ziran_cao_314_gwas.vcf.gz --regions 1:3732143-3735071 > ziram_popuiation_PSY_snp.vcf.gz #1:3732143-3735071表示提取1号染色体这个区间内的变异位点，可以不加.gz

4.利用plink软件将VCF转换为ped和map格式文件

vim vcf_plink.sh#创建一个脚本

# 使用PLINK处理VCF文件的脚本
plink --vcf ziram_popuiation_PSY_snp.vcf \
--maf 0.01 \
--geno 0.2 \
--allow-extra-chr \
--recode \
--out psy_SNP_filtered2

5.map和ped文件整理，使其符合haploview软件输入格式

awk '{print $2 "\t" $4}' psy_SNP_filtered2.map > psy_SNP_filtered2_edited.info #提取map文件中的第二和第四列，形成一个和ped文件名字一样的.info文件。


awk '{gsub(/-9/, "0"); print}' psy_SNP_filtered2.ped > psy_SNP_filtered2_edited.ped #将-9全部替换为0

二，Haploview软件使用

1.安装软件：下载Linux版本的Haploview4.1软件，单独放在一个创建的英文命名的文件夹下

同时需要安装Java环境，安装之后win +r 输入cmd进入命令行，输入java -jar Haploview.jar之后即可以运行Haploview4.1。

2.将整理好的文件输入进去即可

三，单倍型显著性分析

画出单倍型block图后知道有哪些可能存在的单倍型，但haploview软件没有具体的群体单倍型数据，需要做各个单倍型之间某个表型显著性分析时不能进行，需要使用plink进行分析单倍型然后提取出对应基因型数据转化为VCF文件去使用R进行分析作图。

1.plink分析群体中各个株系对应单倍型

plink --file psy_SNP_filtered2 --allow-extra-chr --allow-no-sex --blocks no-pheno-req #no-pheno-req强行指定不要表型，psy_SNP_filtered2是前面分析block的原始ped和map数据。
cat plink.blocks.det #查看有多少block

2.提取出单倍型基因型数据并转为vcf格式文件

提取候选基因中的block区间基因型

plink --file psy_SNP_filtered2 --chr 1 --from-bp 3733853 --to-bp 3733876 --recode --out df8

cat df8.map

#转换为VCF格式
plink --file df7 --recode vcf-iid --out df8_vcf

3.R进行单倍型鉴定及可视化

rm(list = ls())
setwd("D:/研究生期间相关文件/群体基因型及表型数据/CX/genotype/psy相关/")
## 导入基因型vcf数据
library(geneHapR)

vcf8 = import_vcf("df8_vcf.vcf")

# 单倍型分型
hapResult1 <- vcf2hap(vcf8)#默认删除杂合位点
hapResult2 <- vcf2hap(vcf8,hetero_remove = F)#保留杂合位点

# 对单倍型结果进行汇总整理
hapPrefix <- "Hap"        # 单倍型名称的前缀
hapSummary1 <- hap_summary(hapResult1, hapPrefix = hapPrefix)
hapSummary2 <- hap_summary(hapResult2, hapPrefix = hapPrefix)

#保存
write.csv(hapResult1,"df8-1-hapresult.csv")
write.csv(hapResult2,"df8-2-hapresult.csv")
write.csv(hapSummary1,"df8-1-hapSummary.csv")
write.csv(hapSummary2,"df8-2-hapSummary")

#导入表型
phe <- import_AccINFO("phe_ca.txt")

#可视化
geneID = "IbPSY"
#1表格型
p1 <-  plotHapTable(hapSummary1,angle = 45,title = geneID)

p2 <-  plotHapTable(hapSummary2,angle = 45,title = geneID)

print(p1)
print(p2)
# 保存图
ggsave("df7-1-hap.jpg", 
       plot = p1, 
       width = 8, 
       height = 6, 
       dpi = 300)
ggsave("df7-2-hap.jpg", 
       plot = p2, 
       width = 8, 
       height = 6, 
       dpi = 300)

#2不同单倍型间检验
p3 <- hapVsPheno(hap = hapResult1,pheno = phe)
p4 <- hapVsPheno(hap = hapResult2,pheno = phe)

print(p4)
# 保存图
ggsave("df8-test1-hap.jpg", 
       plot = p3$fig_Violin, 
       width = 8, 
       height = 6, 
       dpi = 300)