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 <haplotypes> --map <mapfile> --out <outfile>
我们需要准备 --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 <haplotypes> --map <mapfile> --out <outfile>
现在直接计算就行。
vcf 文件是 phasing 后的压缩文件；
map 文件是每个位点的遗传位置，包括四列：染色体、位点ID、遗传位置、物理位置。
3. 结果文件解读
<locusID > <physicalPos > <’1’ freq > <ihh1 > <ihh0 > < unstandardized iHS >
结果文件包括 6 列，分别是以上内容，而后就按自己的要求来找受选择的区域了。
欢迎评论交流，批评指正！
