Integrated Haplotype Score (iHS) (Voight et al., 2006) 是基于单倍型的选择信号检测软件。iHS 通过比较单倍型频率的衰减速度来检测选择信号。具体来说，iHS 衡量了与一个给定等位基因相关的单倍型的连锁不平衡（LD）在追溯时的变化率。如果一个基因座在群体中受到了强烈的正向选择，那么与该基因座相关的单倍型将显示出不平衡的追溯，即高频等位基因的单倍型将保持相对较高的频率。

GitHub - szpiech/selscan: Haplotype based scans for selection

(https://github.com/szpiech/selscan)

从上面链接可以下载软件 selscan，这个软件可以计算 iHS、EHH、XP-EHH、iHH12、nSL、XP-nSL，本文着重介绍 iHS 的计算过程，其余的看情况也会尝试。

selscan 软件是免编译的，下载可以直接使用。使用方法为：

selscan --ihs --hap --map --out

我们需要准备 --hap 文件和 --map 文件。

通常，我们很难获得单倍型的文件，重测序数据的 VCF 文件也可以，但在此之前需要对文件进行 phasing。

1.1 VCF 文件预处理

基于单倍型的计算只需要种对中一个群体的数据，计算的是种内的单倍型情况.

从总 VCF 文件中调出单个群体的 VCF 文件，XX.list 是该群体的样本列表：

bcftools view XX.vcf.recode.vcf -Oz -o XX.vcf.gz

bcftools index XX.vcf.gz

bcftools view -S   XX.list   XX.vcf.gz > XX.vcf

而后，按染色体将文件分割：

#!/bin/bash

chr_list=("chr01" "chr02" "chr03" "chr04" "chr05" "chr06" "chr07" "chr08" "chr09" "chr10" "chr11" "chr12")

for i in "${chr_list[@]}"

do

    vcftools --vcf nNEP.vcf --chr "$i" --recode --out "$i.vcf"

done

使用 vcftools 对文件进行按染色体分割，每条染色体分开计算。

这部分基本就可以进行 phasing 了，但我在计算的时候，不知道什么原因，部分缺失位点的基因型，理论上应该是  ./.   ，但部分位点变成了   .    ，这会导致 phasing 过程无法读取文件，我也没查到怎么解决，就直接写了个 perl 脚本将  .   的基因型改为了  ./.   。

脚本见 ：修改 VCF 文件中错误编码的基因型

这一段可以忽略，大家可能不会遇到这个问题。

1.2 Beagle 软件进行 phasing

软件官网：

http://faculty.washington.edu/browning/beagle/

同样，下载下来就可以直接使用，不需要编译安装。

java   -Xmx10g   -jar      beagle.18May20.d20.jar      gt=XX.vcf    map=XX.map      out=xx.phased         nthreads=6

上述是计算示例，分别对每条染色体进行计算，得到的结果是一个压缩文件，不需要解压，后面直接读取即可。

gt 指定的是基因型文件，即 VCF 文件；但这里需要指定一个 map 文件，如果是模式植物的话比较好说，给一个遗传图谱即可，物理距离对应遗传距离。

第一列是染色体；第二列是位点ID，任意指定；第三列是遗传位置；第四列是物理位置。

如果没有遗传图谱也可以，默认按 1Mb 为 1cM 进行计算。

1.3 计算位点的遗传距离

在 selscan 计算的时候仍需要一个 map 文件，但这个要求每个位点的遗传距离。

如果没有遗传图谱，可以根据物种总的遗传距离算每个 bp 平均的遗传距离，而后使用物理距离推算。如果有的话，也需要根据遗传图谱中的标记推断每个 SNP 的遗传位置。

标记之间插值填补（interpolate_data.pl）

我写了一个 perl 脚本对遗传图谱中两个标记之间的间隙进行插值填补，获得所有位点的遗传距离。而后可以根据自己的文件中的位点，生成一个 map 文件，进行计算。

perl      interpolate_data.pl    chr01.map    chr01.map.out

zcat   XX_phased.vcf.gz   |   awk '{print $1, $2}' | sed '1,9d' >   chr01 

# 获得VCF文件中所有位点的信息（去除表头）

perl     join.pl     chr01     chr01.map.out     chr01_new.map

# 使用我自己写的一个 perl 脚本确定 VCF 文件中每个位点的遗传位置

awk '{print $1, ".", $3, $2}'   chr01_new.map   >  chr01_new.map.out

# 得到需要的 map 文件格式。

比较两个文件内容，匹配的行写入新的文件（join.pl ）

2. 计算 iHS

selscan --ihs --vcf    --map --out

现在直接计算就行。

vcf 文件是 phasing 后的压缩文件；

map 文件是每个位点的遗传位置，包括四列：染色体、位点ID、遗传位置、物理位置。

3. 结果文件解读

<’1’ freq > < unstandardized iHS >

结果文件包括 6 列，分别是以上内容，而后就按自己的要求来找受选择的区域了。

欢迎评论交流，批评指正！