1、计算等位基因和基因型频率

通过hiblup -h查看帮助可以检索到等位基因频率和基因型频率计算的相关参数--allele-freq和--geno-freq。

1.1、计算等位基因频率

等位基因是指位于一对同源染色体相同位置上控制同一性状不同形态的基因。不同的等位基因将会产生不同的遗传特征。等位基因频率就是某种等位基因占某个种群中全部基因数的比例。使用HIBLUP计算等位基因频率时，命令行输入如下：

hiblup --bfile demo --allele-freq --out demo

--bfile：指定输入的二进制格式的基因型文件的文件名，不包含文件名后缀。hiblup接受的文件格式是plink二进制格式，例如测试数据集中提供的demo.bed、demo.fam、demo.bim文件。

demo.bim：储存每个SNP位点的相关信息，共六列，记录了SNP所在的染色体编号、SNP名称、位于基因组上的位置、碱基对的坐标和等位基因类型。

demo.fam：储存个体家系信息，共六列，分别为家系编号、个体编号、父系编号、母系编号、性别和表型值，其中0和-9表示信息缺失。

demo.bed：以特定的压缩格式储存的基因型信息，在文本编辑器中不可读，详细解释可见：https://www.cog-genomics.org/plink/1.9/formats#bed

如果文件格式为VCF、HapMap、PED/MAP等，可以使用PLINK2.0(https://www.cog-genomics.org/plink/2.0/)、TASSEL(https://www.maizegenetics.net/tassel)、VCFtools(https://vcftools.github.io/man_latest.html#AUTHORS)等软件将格式转换为plink二进制格式后被HIBLUP接受。

①使用TASSEL 将HapMap格式转换成VCF格式：

run_pipeline.pl -SortGenotypeFilePlugin -inputFile demo.hmp.txt -outputFile demo.sort.hmp.txt -fileType Hapmap

run_pipeline.pl -fork1 -h demo.sort.hmp.txt -export-exportType VCF -runfork1

②使用PLINK将VCF格式转换成PLINK二进制格式：

plink2 --vcf demo.vcf --make-bed --out demo

③使用PLINK将PED/MAP格式转换成PLINK二进制格式：

plink2 --ped demo.ped --map demo.map --make-bed --out demo

此外，由于HIBLUP没有填充功能，因此基因型文件中缺失的基因型都将被视为杂合子，在分析中被强制编码为1。为了保证分析的准确性，在使用HIBLUP之前可以使用Beagle(https://faculty.washington.edu/browning/beagle/beagle.html)或其他软件对基因型进行填充。

程序运行效果如下：

本次程序运行计算了demo基因型文件中800个个体1000个SNP位点上的等位基因频率，并把结果储存在demo.afreq文件中，程序运行的信息记录在demo.log文件中。

查看demo.afreq文件，格式如下：

第一列为SNP位点的名称；第二列和第三列为该SNP位点上的两个等位基因a1、a2的形式；第四列为a1等位基因的频率。a2等位基因的频率可以通过1-a1等位基因的频率得到。

1.2、计算基因型频率

基因型是指等位基因在某个基因座上的不同组合形式。基因型频率就是群体中某一种特定的基因型个体数在群体中所有个体所占的比例。使用HIBLUP计算基因型频率时，命令行输入如下：

hiblup --bfile demo --geno-freq --out demo

--bfile：指定输入的二进制格式的基因型文件的文件名，不包含文件名后缀，即demo.bed、demo.fam、demo.bim；

--geno-freq：计算基因型频率；

--out：指定输出的文件路径和文件名。

程序运行效果如下：

这里计算了demo基因型文件中800个个体1000个SNP位点的基因型频率，结果储存在demo.gfreq文件中，程序运行的信息记录在demo.log文件中。

查看demo.gfreq文件，格式如下：

第一列为SNP位点的名称；第二列和第三列分别为该SNP位点上的两个等位基因a1、a2的形式；第四列和第五列分别为a1a1基因型和a2a2基因型的频率。a1a2基因型的频率可以通过1-a1a1基因型频率-a2a2基因型频率得到。

HIBLUP还提供了进行SNP位点筛选的参数--extract和--exclude。所需要输入的文件格式如snp.filter.txt所示，为一行一个的SNP位点名称。

--extract：输入包含分析中需要使用的SNP名称的文件，分析时将只对文件中包含的SNP位点进行计算；

--exclude：输入包含分析中需要删除的SNP名称的文件，分析时将删除文件中包含的SNP位点，对其他位点进行计算。

2、计算个体纯合率和杂合率

HIBLUP提供了用于计算个体纯合率和杂合率的参数--homo和—hete。

2.1、计算个体纯合率

纯合子是同源染色体上相同位点等位基因相同的基因型。在个体中，纯合率是指由两个相同等位基因组成的纯合基因座占该个体所有基因座中的比例。使用HIBLUP计算纯合率时，命令行输入如下：

hiblup --bfile demo --homo --out demo

--bfile：指定输入的二进制格式的基因型文件的文件名，不包含文件名后缀；

--homo：计算个体纯合率；

--out：指定输出的文件路径和文件名。

程序运行效果如下：

运行完成后生成demo.homo文件和demo.log文件。基因组文件中800个个体在1000个基因座中a1或a2等位基因纯合的频率保存在demo.homo文件中。

查看demo.homo文件，格式如下：

第一列为个体ID；第二列和第三列分别为该个体a1等位基因纯合和a2等位基因纯合的频率。两列相加可以得到个体总纯合率。杂合率可以通过1-a1等位基因纯合频率-a2等位基因纯合频率得到，也可以通过--hete参数直接计算，计算方法如下。

2.2、计算个体杂合率

杂合子是同源染色体同一位点上的两个等位基因不相同的基因型。在个体中，杂合率是指由两个不同的等位基因组成杂合基因座占该个体所有基因座中的比例。使用HIBLUP计算杂合率时，命令行输入如下：

hiblup --bfile demo --hete --out demo

--bfile：指定输入的二进制格式的基因型文件的文件名，不包含文件名后缀；

--hete：计算个体杂合率；

--out：指定输出的文件路径和文件名。

程序运行效果如下：

基因组文件中800个个体在1000个基因座中a1a2等位基因杂合的频率保存在demo.hete文件中。

查看demo.hete文件，格式如下：

第一列为基因型个体ID；第二列为该个体a1a2等位基因杂合的频率。

HIBLUP还提供了进行个体筛选的参数--keep和--remove。所需要输入的文件格式如id.filter.txt所示，为一行一个的个体名称。

--keep：输入包含分析中需要使用的个体ID的文件，分析时将只对文件中包含的个体进行计算；

--remove：输入包含分析中需要删除的个体ID的文件，分析时将删除文件中包含的个体，对其他个体进行计算。

3、修改线程数

为了达到更高的计算效率，提高资源利用率，HIBLUP会默认获取OpenMP环境变量中最大线程数用于分析。由于没有指定线程数，这里默认使用32个线程进行计算。

如果需要指定线程数，只需要在计算时使用--threads参数指定用于分析的线程数。比如将上述计算基因型频率的命令行代码更换为：

hiblup --bfile demo --geno-freq –threads 10 --out demo

就可以只使用10个线程进行计算，但计算效率也可能会受到影响。