大家好，我是鄧飛，今天介紹一下Haploview分析的方法。

本文詳細(xì)介紹了Haploview單倍型的分析方法，從背景介紹、軟件安裝、數(shù)據(jù)準(zhǔn)備、作圖分析以及批量操作，通過(guò)學(xué)習(xí)和實(shí)踐操作，讀者可以輕松掌握單倍型分析的方法。

一、單倍型介紹

為何要做單倍型分析?

我們做完GWAS分析，得到了顯著性位點(diǎn)，注釋到了上下游的基因，這時(shí)，一個(gè)想法浮現(xiàn)在眼前：你如何證明你找到的基因不是假陽(yáng)性？？？

答案就是單倍型分析，看一下顯著性位點(diǎn)附近的區(qū)域，是否處于一個(gè)高度連鎖的區(qū)域（block），看一下基因是否在block里面，如果顯著位點(diǎn)附近有高連鎖的BLOCK并且注釋的基因也在block里面，可以證明挖掘的基因沒(méi)問(wèn)題，結(jié)果八九不離十了，十分可靠。 

那如何做單倍型分析呢？  

如果按照分析思路的話，是選擇顯著性為點(diǎn)上下游的區(qū)域，計(jì)算SNP之間的LD值，然后根據(jù)某個(gè)閾值進(jìn)行劃分Block，如果有block，那么block區(qū)域內(nèi)只有少數(shù)的組合，這些少數(shù)的組合就是單倍型。我們定位基因，或者分子標(biāo)記輔助，都會(huì)用到單倍型。 

好消息是，不用自己手動(dòng)計(jì)算LD值，然后變成劃分block了，有現(xiàn)成的軟件。壞消息是軟件也要學(xué)習(xí)，目前主流的兩款軟：Haploview和LDblockshow，前者是桌面版軟件，后者是命令行軟件，兩者結(jié)果基本一致。

二、Haploview軟件安裝

1. 軟件下載

「官網(wǎng)：」 https://www./haploview/downloads#JAR

「windows系統(tǒng)：」

「Linux系統(tǒng)：」

2. 配置java環(huán)境

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

下載安裝好之后，在終端運(yùn)行java，出現(xiàn)幫助文檔，說(shuō)明配置成功。

3 windows系統(tǒng)

cmd終端打開(kāi)，鍵入java，出現(xiàn)下面界面，說(shuō)明配置成功。

找到快捷方式：

打開(kāi)軟件：

4 Linux系統(tǒng)

終端打開(kāi)，鍵入java，出現(xiàn)下面界面，說(shuō)明配置成功。

終端運(yùn)行，即可出現(xiàn)圖形界面。

注意：Linux系統(tǒng)，需要窗口界面，如果是終端登錄，可以安裝xming，顯示窗口界面。

5 windows報(bào)錯(cuò)解決方案

注意，有時(shí)候會(huì)有問(wèn)題，windows系統(tǒng)安裝Haploview總是錯(cuò)誤，解決方案是直接用jar文件，在終端cmd中運(yùn)行：

有個(gè)小伙伴說(shuō)，Haploview軟件在windows系統(tǒng)中安裝報(bào)錯(cuò)，說(shuō)是配置了java和javac環(huán)境，還是不錯(cuò)，我是不信的。 

于是我用自己的業(yè)余電腦試了一下，果然報(bào)錯(cuò)：

安裝報(bào)錯(cuò)信息：

配置好的java  

然后導(dǎo)入java.exe文件： 

太難了。重頭試了一遍，還是報(bào)錯(cuò)。

想到一個(gè)方法，既然Linux和Mac可以運(yùn)行Haploview.jar文件，那么windows系統(tǒng)也可以，測(cè)試一下： 

然后就可以運(yùn)行了。 

寫(xiě)一個(gè)bat腳本： 

雙擊一下，就有Haploview的GUI界面了。

三、單倍型數(shù)據(jù)準(zhǔn)備

1. 數(shù)據(jù)準(zhǔn)備

需要做單倍型分析的是基因型數(shù)據(jù)，一般是顯著性的SNP，提取上下游500kb，然后進(jìn)行block的分析。

這里，準(zhǔn)備的是plink數(shù)據(jù)，比如我們要提?。?/p>

• 染色體是6

• 開(kāi)始位置是1000000

• 終止位置是2000000

將其轉(zhuǎn)化為plink的map和ped數(shù)據(jù)：

2. 整理數(shù)據(jù)

將map的第二列和第四列提取出來(lái)，保存為a1.info文件。

ped數(shù)據(jù)，保持不變：

3. 導(dǎo)入數(shù)據(jù)

選擇第一種格式：Linkage Format，然后將ped數(shù)據(jù)導(dǎo)入到Data File中，將info數(shù)據(jù)導(dǎo)入到Locus Information File文件中。

結(jié)果：

查看Block：

查看TaggerSNP：

上面就是下數(shù)據(jù)分析實(shí)操方法。

4. plink格式整理常見(jiàn)錯(cuò)誤

有些朋友用自己的數(shù)據(jù)分析時(shí)，會(huì)有各種問(wèn)題，最近星球上有小伙伴發(fā)了一個(gè)帖子，敘述了自己的問(wèn)題，各種嘗試，還是錯(cuò)誤，淡淡的憂傷和砸電腦的沖動(dòng)……  

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

作為分析數(shù)據(jù)的老司機(jī)，新手遇到的錯(cuò)誤都會(huì)遇到過(guò)，而且也有解決方案，這里總結(jié)一下常見(jiàn)的錯(cuò)誤。 

1，vcf變?yōu)閜link，或者excel變?yōu)閜link格式，然后直接讀入到Haploview，卒。

vcf變?yōu)閜link，以及plink提取snp導(dǎo)出plink格式，或者excel的格式變?yōu)閜link格式，這些數(shù)據(jù)都建議用plink重新跑一下，確保數(shù)據(jù)沒(méi)問(wèn)題。   如果plink運(yùn)行報(bào)錯(cuò)，就不要往下走了，先把這個(gè)問(wèn)題解決掉！ 

2，plink格式的map數(shù)據(jù)，需要變?yōu)閕nfo格式，簡(jiǎn)單來(lái)說(shuō)，就是提取map的第二列和第四列，生成info為后綴的文件  

3，ped中的表型數(shù)據(jù)默認(rèn)為-9，需要改為0。這個(gè)是必須的，否則就報(bào)錯(cuò)。 

4，另外，還有一個(gè)坑，ped的ID編號(hào)中不能有-

這里建議用下劃線代替，搞定這些后，就可以用Haploview讀取數(shù)據(jù)，進(jìn)行單倍型分析了。

代碼匯總：

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

還有一個(gè)靈魂拷問(wèn)的問(wèn)題，不要用所有的基因型數(shù)據(jù)做單倍型分析，只需要把顯著SNP上下游一段距離（比如50kb）的位點(diǎn)提取出來(lái)，進(jìn)行單倍型分析。

四、Haploview結(jié)果解讀

1. LDblock整體結(jié)果

上圖就是最常見(jiàn)的LDblock，該圖的結(jié)果解讀。

2. SNP在染色體的分布

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

3. SNP的名稱信息

中間是SNP的名稱，用細(xì)線聯(lián)系在一起

4. block解釋

最下面紅白的正方形是LD值的可視化，每一個(gè)正方形是兩兩SNP的LD結(jié)果，顏色越淡說(shuō)明LD值越小，如果相鄰的SNP之間的LD大于某個(gè)閾值（比如0.9），那么就構(gòu)成一個(gè)block，下圖中的兩個(gè)紅框里面的黑框，就是兩個(gè)LDblock，第一個(gè)block包括的SNP有10,11,12三個(gè)SNP，block的距離為82kb，第二個(gè)block包括兩個(gè)snp，包括14和15兩個(gè)snp，block的距離為32kb。

5. 查看block的頻率和他們之間的聯(lián)系

下圖中，第一個(gè)block中，一共三個(gè)SNP，單倍型分別是：TTC，TTA，CCA，TCA，他們的頻率分別是0.548,0.281,0.09和0.078，它們的頻率之和為1。第二個(gè)block一共有兩個(gè)SNP，單倍型分別是AG，GA和AA，頻率分別是0.402,0.565,0.034，他們之間的頻率之和為1.

最下面的0.67是兩個(gè)block的關(guān)聯(lián)，兩個(gè)block的線是兩者的關(guān)聯(lián)性，線條越黑，說(shuō)明關(guān)聯(lián)性越強(qiáng)。

6. 查看TaggerSNP

這里有兩個(gè)block，可以選擇兩個(gè)TaggerSNP代表這兩個(gè)block

五、LDBlockshow軟件安裝

「要實(shí)現(xiàn)的下面的圖：」

• 最下方的熱圖是兩兩SNP之間的LD值，越高越紅，比較紅的區(qū)域構(gòu)成一個(gè)Block（用黑線連起來(lái)）

• 如果提供gff文件，可以顯示基因的上游、下游、外顯子、內(nèi)含子區(qū)域

• 上面是位點(diǎn)的曼哈頓圖，是區(qū)域性的曼哈頓圖

• 位點(diǎn)之間，也可以根據(jù)LD值進(jìn)行可視化，以最顯著的位點(diǎn)為四方形，其它位點(diǎn)與其LD值的大小呈現(xiàn)不同的顏色

軟件介紹：github鏈接：https://github.com/BGI-shenzhen/

• A:整體上宏觀上用：

  PopLDdecay 軟件 ，軟件己經(jīng)生物信息Bioinformatics雜志發(fā)表online（見(jiàn)我的學(xué)習(xí)筆記：

  LD衰減圖繪制--PopLDdecay）

• B: 從局部上查看用：

  LDBlockShow軟件， 軟件已經(jīng)正式被 briefings in bioinformatics （影響分子8.99）的雜志接收

1. 數(shù)據(jù)準(zhǔn)備

• vcf格式的數(shù)據(jù)，InVCF

• plink二進(jìn)制文件，InPlink

• plink文本文件，InPlink

2. 軟件安裝

網(wǎng)址：https://github.com/BGI-shenzhen/LDBlockShow

中文說(shuō)明書(shū)：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. 軟件測(cè)試

數(shù)據(jù)：file.vcf

代碼：這里，繪制染色體1，位置區(qū)間是：49670000:49780000

結(jié)果：

4. 進(jìn)階：Heatmap + block

vcf文件：Test.vcf.gz

命令：

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

結(jié)果文件：

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

5. 進(jìn)階：Heatmap + block + GWAS

考慮GWAS的結(jié)果，加入?yún)?shù)：-InGWAS gwas.pvalue

vcf文件：Test.vcf.gz

GWAS結(jié)果文件：三列，Chr, Position, Pvalue，沒(méi)有行頭

$ 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-InVCFTest.vcf.gz-OutPutre2-Regionchr11:24100000:24200000-InGWASgwas.pvalue-OutPng-SeleVar 1

結(jié)果：

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

結(jié)果中包括熱圖，block圖和GWAS圖合并起來(lái)了。

上面的圖，可以通過(guò)ShowLDSVG軟件，進(jìn)一步優(yōu)化：

• -Cutline，閾值定義為7

• -ShowNum，顯示LD值

• -PointSize，顯示點(diǎn)大小

結(jié)果：

6. Heatmap + block + GWAS + Annotation

相比較上圖，增加了注釋的信息。

文件需要：

• vcf，vcf格式的文件    

• gwas_pvalue，三列的gwas結(jié)果（Chr，Position，Pvalue），無(wú)行頭

• gff文件，注釋文件

命令：

LDBlockShow-InVCFTest.vcf.gz-OutPutre3-Regionchr11:24100000:24200000-InGWASgwas.pvalue-OutPng-SeleVar 1 -InGFFIn.gff

也可以增加SNP的名稱：

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

命令：

7. 進(jìn)階：LDblock+GWAS+Annotation+Locuszoom

可以通過(guò)-TopSite在GWAS圖中顯示最顯著位點(diǎn)與其它位點(diǎn)的LD關(guān)系。

下圖中，最顯著的位點(diǎn)為四邊形，其它顏色，紅色表示LD高，其它顏色表示LD低。在上圖的基礎(chǔ)上，增加了最顯著位點(diǎn)與其它位點(diǎn)的LD情況。

六、批量繪制單倍型圖

有小伙伴在星球詢問(wèn)GWAS分析的顯著性位點(diǎn)，如何批量繪制LDblock，之前寫(xiě)過(guò)LDblock的安裝和使用說(shuō)明：（LDblock繪制連鎖不平衡和單體型圖）

問(wèn)題有3個(gè)核心點(diǎn)： 

1，snp文件總是有^M符號(hào)，怎么去除
2，如何批量對(duì)顯著位點(diǎn)進(jìn)行LDblock分析
3，有沒(méi)有現(xiàn)成代碼

今天就搞定這個(gè)問(wèn)題。 首先，看一下正常LDblock分析的代碼： 

有vcf文件，染色體，物理位置區(qū)間，結(jié)果文件名稱。 

來(lái)一個(gè)腳本： 

$ cat ~/bin/plot_ld_plot_vcf.sh #!/bin/bash
if [ $# != 4 ]then    echo $0 name.vcf chronome_number location 20000; 第一個(gè)是vcf或者vcf.gz文件，第二個(gè)是染色體位置，第三個(gè)是物理位置，第四個(gè)是上下游區(qū)間    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文件，染色體，物理位置，上下游區(qū)間，就會(huì)自動(dòng)分析。 

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

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

上下游區(qū)間是100kb，那么先將每個(gè)snp生成一個(gè)命令行，然后運(yùn)行就行了。 

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

然后運(yùn)行temp.sh文件就可以了。靜等結(jié)果就可以了。