从入门到出家：单倍型Haploview分析（万字详解）

大家好，我是邓飞，今天介绍一下Haploview分析的方法。

本文详细介绍了Haploview单倍型的分析方法，从背景介绍、软件安装、数据准备、作图分析以及批量操作，通过学习和实践操作，读者可以轻松掌握单倍型分析的方法。

一、单倍型介绍

为何要做单倍型分析?

我们做完GWAS分析，得到了显著性位点，注释到了上下游的基因，这时，一个想法浮现在眼前：你如何证明你找到的基因不是假阳性？？？

答案就是单倍型分析，看一下显著性位点附近的区域，是否处于一个高度连锁的区域（block），看一下基因是否在block里面，如果显著位点附近有高连锁的BLOCK并且注释的基因也在block里面，可以证明挖掘的基因没问题，结果八九不离十了，十分可靠。 

那如何做单倍型分析呢？  

如果按照分析思路的话，是选择显著性为点上下游的区域，计算SNP之间的LD值，然后根据某个阈值进行划分Block，如果有block，那么block区域内只有少数的组合，这些少数的组合就是单倍型。我们定位基因，或者分子标记辅助，都会用到单倍型。 

好消息是，不用自己手动计算LD值，然后变成划分block了，有现成的软件。坏消息是软件也要学习，目前主流的两款软：Haploview和LDblockshow，前者是桌面版软件，后者是命令行软件，两者结果基本一致。

二、Haploview软件安装

1. 软件下载

「官网：」 https://www.broadinstitute.org/haploview/downloads#JAR

「windows系统：」

「Linux系统：」

2. 配置java环境

https://www.java.com/zh-CN/download/

下载安装好之后，在终端运行java，出现帮助文档，说明配置成功。

3 windows系统

cmd终端打开，键入java，出现下面界面，说明配置成功。

找到快捷方式：

打开软件：

4 Linux系统

终端打开，键入java，出现下面界面，说明配置成功。

终端运行，即可出现图形界面。

注意：Linux系统，需要窗口界面，如果是终端登录，可以安装xming，显示窗口界面。

5 windows报错解决方案

注意，有时候会有问题，windows系统安装Haploview总是错误，解决方案是直接用jar文件，在终端cmd中运行：

有个小伙伴说，Haploview软件在windows系统中安装报错，说是配置了java和javac环境，还是不错，我是不信的。 

于是我用自己的业余电脑试了一下，果然报错：

安装报错信息：

配置好的java  

然后导入java.exe文件： 

太难了。重头试了一遍，还是报错。

想到一个方法，既然Linux和Mac可以运行Haploview.jar文件，那么windows系统也可以，测试一下： 

然后就可以运行了。 

写一个bat脚本： 

双击一下，就有Haploview的GUI界面了。

三、单倍型数据准备

1. 数据准备

需要做单倍型分析的是基因型数据，一般是显著性的SNP，提取上下游500kb，然后进行block的分析。

这里，准备的是plink数据，比如我们要提取：

• 染色体是6

• 开始位置是1000000

• 终止位置是2000000

将其转化为plink的map和ped数据：

2. 整理数据

将map的第二列和第四列提取出来，保存为a1.info文件。

ped数据，保持不变：

3. 导入数据

选择第一种格式：Linkage Format，然后将ped数据导入到Data File中，将info数据导入到Locus Information File文件中。

结果：

查看Block：

查看TaggerSNP：

上面就是下数据分析实操方法。

4. plink格式整理常见错误

有些朋友用自己的数据分析时，会有各种问题，最近星球上有小伙伴发了一个帖子，叙述了自己的问题，各种尝试，还是错误，淡淡的忧伤和砸电脑的冲动……  

https://wx.zsxq.com/dweb2/index/group/48844515845858

作为分析数据的老司机，新手遇到的错误都会遇到过，而且也有解决方案，这里总结一下常见的错误。 

1，vcf变为plink，或者excel变为plink格式，然后直接读入到Haploview，卒。

vcf变为plink，以及plink提取snp导出plink格式，或者excel的格式变为plink格式，这些数据都建议用plink重新跑一下，确保数据没问题。   如果plink运行报错，就不要往下走了，先把这个问题解决掉！ 

2，plink格式的map数据，需要变为info格式，简单来说，就是提取map的第二列和第四列，生成info为后缀的文件  

3，ped中的表型数据默认为-9，需要改为0。这个是必须的，否则就报错。 

4，另外，还有一个坑，ped的ID编号中不能有-

这里建议用下划线代替，搞定这些后，就可以用Haploview读取数据，进行单倍型分析了。

代码汇总：

plink --file file --maf 0.01 --geno 0.3 --recode --out qc300sed  -i 's/-/_/g' qc300.pedsed  -i 's/-9/0/g' qc300.pedawk '{print $2,$4}' qc300.map >qc300.info

还有一个灵魂拷问的问题，不要用所有的基因型数据做单倍型分析，只需要把显著SNP上下游一段距离（比如50kb）的位点提取出来，进行单倍型分析。

四、Haploview结果解读

1. LDblock整体结果

上图就是最常见的LDblock，该图的结果解读。

2. SNP在染色体的分布

最上面是SNP的物理位置，有些是均匀的，有些是不均匀的

3. SNP的名称信息

中间是SNP的名称，用细线联系在一起

4. block解释

最下面红白的正方形是LD值的可视化，每一个正方形是两两SNP的LD结果，颜色越淡说明LD值越小，如果相邻的SNP之间的LD大于某个阈值（比如0.9），那么就构成一个block，下图中的两个红框里面的黑框，就是两个LDblock，第一个block包括的SNP有10,11,12三个SNP，block的距离为82kb，第二个block包括两个snp，包括14和15两个snp，block的距离为32kb。

5. 查看block的频率和他们之间的联系

下图中，第一个block中，一共三个SNP，单倍型分别是：TTC，TTA，CCA，TCA，他们的频率分别是0.548,0.281,0.09和0.078，它们的频率之和为1。第二个block一共有两个SNP，单倍型分别是AG，GA和AA，频率分别是0.402,0.565,0.034，他们之间的频率之和为1.

最下面的0.67是两个block的关联，两个block的线是两者的关联性，线条越黑，说明关联性越强。

6. 查看TaggerSNP

这里有两个block，可以选择两个TaggerSNP代表这两个block

五、LDBlockshow软件安装

「要实现的下面的图：」

• 最下方的热图是两两SNP之间的LD值，越高越红，比较红的区域构成一个Block（用黑线连起来）

• 如果提供gff文件，可以显示基因的上游、下游、外显子、内含子区域

• 上面是位点的曼哈顿图，是区域性的曼哈顿图

• 位点之间，也可以根据LD值进行可视化，以最显著的位点为四方形，其它位点与其LD值的大小呈现不同的颜色

软件介绍：github链接：https://github.com/BGI-shenzhen/

• A:整体上宏观上用：

  PopLDdecay 软件 ，软件己经生物信息Bioinformatics杂志发表online（见我的学习笔记：

  LD衰减图绘制--PopLDdecay）

• B: 从局部上查看用：

  LDBlockShow软件， 软件已经正式被 briefings in bioinformatics （影响分子8.99）的杂志接收

1. 数据准备

• vcf格式的数据，InVCF

• plink二进制文件，InPlink

• plink文本文件，InPlink

2. 软件安装

网址：https://github.com/BGI-shenzhen/LDBlockShow

中文说明书：https://github.com/hewm2008/LDBlockShow/blob/main/LDBlockShow_Manual_Chinese.pdf

安装代码：

git clone https://github.com/hewm2008/LDBlockShow.git   cd LDBlockShow ; chmod 755 configure  ; ./configure;           make;          

3. 软件测试

数据：file.vcf

代码：这里，绘制染色体1，位置区间是：49670000:49780000

结果：

4. 进阶：Heatmap + block

vcf文件：Test.vcf.gz

命令：

LDBlockShow -InVCF Test.vcf.gz -OutPut re1 -Region chr11:24100000:24200000 -OutPng -SeleVar 1

结果文件：

re1.blocks.gz  re1.png  re1.site.gz  re1.svg  re1.TriangleV.gz

5. 进阶：Heatmap + block + GWAS

考虑GWAS的结果，加入参数：-InGWAS gwas.pvalue

vcf文件：Test.vcf.gz

GWAS结果文件：三列，Chr, Position, Pvalue，没有行头

$ head gwas.pvalue   chr11 24142640 0.00009   chr11 24142660 1.02e-9   chr11 24142669 1e-9   chr11 24142692 0.5   chr11 24142724 0.6   chr11 24142756 0.001

命令：

LDBlockShow -InVCF Test.vcf.gz -OutPut re2 -Region chr11:24100000:24200000 -InGWAS gwas.pvalue -OutPng -SeleVar 1

结果：

re2.blocks.gz  re2.png  re2.site.gz  re2.svg  re2.TriangleV.gz

结果中包括热图，block图和GWAS图合并起来了。

上面的图，可以通过ShowLDSVG软件，进一步优化：

• -Cutline，阈值定义为7

• -ShowNum，显示LD值

• -PointSize，显示点大小

结果：

6. Heatmap + block + GWAS + Annotation

相比较上图，增加了注释的信息。

文件需要：

• vcf，vcf格式的文件    

• gwas_pvalue，三列的gwas结果（Chr，Position，Pvalue），无行头

• gff文件，注释文件

命令：

LDBlockShow -InVCF Test.vcf.gz -OutPut re3 -Region chr11:24100000:24200000 -InGWAS gwas.pvalue -OutPng -SeleVar 1 -InGFF In.gff

也可以增加SNP的名称：

$ cat Spe.snp chr11 24142660 chr11 24142669 SpeA chr11 24142760 SpeB

命令：

7. 进阶：LDblock+GWAS+Annotation+Locuszoom

可以通过-TopSite在GWAS图中显示最显著位点与其它位点的LD关系。

下图中，最显著的位点为四边形，其它颜色，红色表示LD高，其它颜色表示LD低。在上图的基础上，增加了最显著位点与其它位点的LD情况。

六、批量绘制单倍型图

有小伙伴在星球询问GWAS分析的显著性位点，如何批量绘制LDblock，之前写过LDblock的安装和使用说明：（LDblock绘制连锁不平衡和单体型图）

问题有3个核心点： 

1，snp文件总是有^M符号，怎么去除
2，如何批量对显著位点进行LDblock分析
3，有没有现成代码

今天就搞定这个问题。 首先，看一下正常LDblock分析的代码： 

有vcf文件，染色体，物理位置区间，结果文件名称。 

来一个脚本： 

$ cat ~/bin/plot_ld_plot_vcf.sh #!/bin/bash
if [ $# != 4 ]then    echo $0 name.vcf chronome_number location 20000; 第一个是vcf或者vcf.gz文件，第二个是染色体位置，第三个是物理位置，第四个是上下游区间    exit 1fi
vcf=$1;chr=$2;pos=$3;qujian=$4st=$(( $pos - $qujian))en=$(( $pos + qujian))
ot=re_${chr}_${pos}echo $chr $pos $pos >temp.txt~/bin/LDBlockShow -InVCF $vcf -Region ${chr}:${st}:${en} -SeleVar 1 -BlockType 2 -NoShowLDist 10000000 -OutPut ${ot}

上面的脚本，给定vcf文件，染色体，物理位置，上下游区间，就会自动分析。 

如何批量处理呢？ 比如snp.txt文件，如下图所示：

1 2807198827 1435474132 1419428142 1419427351 17550492610 336195821 24513624410 1340346861 540492201 2852738455 182405365 832334878 1232729356 9326360510 14417231610 35982628 47665935 2008132768 160536300

上下游区间是100kb，那么先将每个snp生成一个命令行，然后运行就行了。 

awk走起： 

awk '{print "bash plot_ld_block_vcf.sh ../geno/genotype.vcf", $1,$2,100000}' snp.txt >temp.sh

生成的脚本文件：

$ cat temp.sh bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 280719882 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 7 143547413 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 2 141942814 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 2 141942735 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 175504926 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 33619582 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 245136244 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 134034686 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 54049220 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 285273845 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 5 18240536 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 5 83233487 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 8 123272935 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 6 93263605 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 144172316 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 3598262 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 8 4766593 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 5 200813276 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 8 160536300 100000

然后运行temp.sh文件就可以了。静等结果就可以了。

想要更好的学习和交流，快来加入飞哥的知识星球，这是一个生物统计+数量遗传学+GWAS+GS的社区，在这里你可以向飞哥提问、帮你制定学习计划、跟着飞哥一起做实战项目，冲冲冲。点击这里加入吧：飞哥的学习圈子

资源推荐：

1，快来领取 | 飞哥的GWAS分析教程

2，飞哥汇总 | 入门数据分析资源推荐

3，数量遗传学，分享几本书的电子版

4，R语言学习看最新版的电子书不香嘛？

5，书籍及配套代码领取--统计遗传分析导论

6，飞哥的学习圈子