1、計(jì)算等位基因和基因型頻率

通過hiblup -h查看幫助可以檢索到等位基因頻率和基因型頻率計(jì)算的相關(guān)參數(shù)--allele-freq和--geno-freq懂版。

1.1鹃栽、計(jì)算等位基因頻率

等位基因是指位于一對(duì)同源染色體相同位置上控制同一性狀不同形態(tài)的基因。不同的等位基因?qū)?huì)產(chǎn)生不同的遺傳特征躯畴。等位基因頻率就是某種等位基因占某個(gè)種群中全部基因數(shù)的比例民鼓。使用HIBLUP計(jì)算等位基因頻率時(shí)，命令行輸入如下：

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

--bfile：指定輸入的二進(jìn)制格式的基因型文件的文件名蓬抄，不包含文件名后綴丰嘉。hiblup接受的文件格式是plink二進(jìn)制格式，例如測(cè)試數(shù)據(jù)集中提供的demo.bed嚷缭、demo.fam饮亏、demo.bim文件。

demo.bim：儲(chǔ)存每個(gè)SNP位點(diǎn)的相關(guān)信息阅爽，共六列路幸，記錄了SNP所在的染色體編號(hào)、SNP名稱优床、位于基因組上的位置劝赔、堿基對(duì)的坐標(biāo)和等位基因類型。

demo.fam：儲(chǔ)存?zhèn)€體家系信息胆敞，共六列着帽，分別為家系編號(hào)、個(gè)體編號(hào)移层、父系編號(hào)仍翰、母系編號(hào)、性別和表型值观话，其中0和-9表示信息缺失予借。

demo.bed：以特定的壓縮格式儲(chǔ)存的基因型信息，在文本編輯器中不可讀，詳細(xì)解釋可見：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)等軟件將格式轉(zhuǎn)換為plink二進(jìn)制格式后被HIBLUP接受。

①使用TASSEL 將HapMap格式轉(zhuǎn)換成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格式轉(zhuǎn)換成PLINK二進(jìn)制格式：

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

③使用PLINK將PED/MAP格式轉(zhuǎn)換成PLINK二進(jìn)制格式：

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

此外狞换，由于HIBLUP沒有填充功能避咆，因此基因型文件中缺失的基因型都將被視為雜合子，在分析中被強(qiáng)制編碼為1修噪。為了保證分析的準(zhǔn)確性查库，在使用HIBLUP之前可以使用Beagle(https://faculty.washington.edu/browning/beagle/beagle.html)或其他軟件對(duì)基因型進(jìn)行填充。

程序運(yùn)行效果如下：

本次程序運(yùn)行計(jì)算了demo基因型文件中800個(gè)個(gè)體1000個(gè)SNP位點(diǎn)上的等位基因頻率黄琼，并把結(jié)果儲(chǔ)存在demo.afreq文件中樊销，程序運(yùn)行的信息記錄在demo.log文件中。

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

第一列為SNP位點(diǎn)的名稱现柠；第二列和第三列為該SNP位點(diǎn)上的兩個(gè)等位基因a1院领、a2的形式弛矛；第四列為a1等位基因的頻率。a2等位基因的頻率可以通過1-a1等位基因的頻率得到比然。

1.2丈氓、計(jì)算基因型頻率

基因型是指等位基因在某個(gè)基因座上的不同組合形式∏糠ǎ基因型頻率就是群體中某一種特定的基因型個(gè)體數(shù)在群體中所有個(gè)體所占的比例万俗。使用HIBLUP計(jì)算基因型頻率時(shí)，命令行輸入如下：

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

--bfile：指定輸入的二進(jìn)制格式的基因型文件的文件名饮怯，不包含文件名后綴闰歪，即demo.bed、demo.fam蓖墅、demo.bim库倘；

--geno-freq：計(jì)算基因型頻率；

--out：指定輸出的文件路徑和文件名论矾。

程序運(yùn)行效果如下：

這里計(jì)算了demo基因型文件中800個(gè)個(gè)體1000個(gè)SNP位點(diǎn)的基因型頻率教翩，結(jié)果儲(chǔ)存在demo.gfreq文件中，程序運(yùn)行的信息記錄在demo.log文件中贪壳。

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

第一列為SNP位點(diǎn)的名稱；第二列和第三列分別為該SNP位點(diǎn)上的兩個(gè)等位基因a1、a2的形式彪笼；第四列和第五列分別為a1a1基因型和a2a2基因型的頻率钻注。a1a2基因型的頻率可以通過1-a1a1基因型頻率-a2a2基因型頻率得到。

HIBLUP還提供了進(jìn)行SNP位點(diǎn)篩選的參數(shù)--extract和--exclude配猫。所需要輸入的文件格式如snp.filter.txt所示队寇，為一行一個(gè)的SNP位點(diǎn)名稱。

--extract：輸入包含分析中需要使用的SNP名稱的文件章姓，分析時(shí)將只對(duì)文件中包含的SNP位點(diǎn)進(jìn)行計(jì)算佳遣；

--exclude：輸入包含分析中需要?jiǎng)h除的SNP名稱的文件，分析時(shí)將刪除文件中包含的SNP位點(diǎn)凡伊，對(duì)其他位點(diǎn)進(jìn)行計(jì)算零渐。

2、計(jì)算個(gè)體純合率和雜合率

HIBLUP提供了用于計(jì)算個(gè)體純合率和雜合率的參數(shù)--homo和—hete系忙。

2.1诵盼、計(jì)算個(gè)體純合率

純合子是同源染色體上相同位點(diǎn)等位基因相同的基因型。在個(gè)體中银还，純合率是指由兩個(gè)相同等位基因組成的純合基因座占該個(gè)體所有基因座中的比例风宁。使用HIBLUP計(jì)算純合率時(shí)，命令行輸入如下：

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

--bfile：指定輸入的二進(jìn)制格式的基因型文件的文件名蛹疯，不包含文件名后綴戒财；

--homo：計(jì)算個(gè)體純合率；

--out：指定輸出的文件路徑和文件名捺弦。

程序運(yùn)行效果如下：

運(yùn)行完成后生成demo.homo文件和demo.log文件饮寞。基因組文件中800個(gè)個(gè)體在1000個(gè)基因座中a1或a2等位基因純合的頻率保存在demo.homo文件中列吼。

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

第一列為個(gè)體ID；第二列和第三列分別為該個(gè)體a1等位基因純合和a2等位基因純合的頻率寞钥。兩列相加可以得到個(gè)體總純合率慌申。雜合率可以通過1-a1等位基因純合頻率-a2等位基因純合頻率得到，也可以通過--hete參數(shù)直接計(jì)算理郑，計(jì)算方法如下蹄溉。

2.2、計(jì)算個(gè)體雜合率

雜合子是同源染色體同一位點(diǎn)上的兩個(gè)等位基因不相同的基因型香浩。在個(gè)體中类缤，雜合率是指由兩個(gè)不同的等位基因組成雜合基因座占該個(gè)體所有基因座中的比例。使用HIBLUP計(jì)算雜合率時(shí)邻吭，命令行輸入如下：

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

--bfile：指定輸入的二進(jìn)制格式的基因型文件的文件名餐弱，不包含文件名后綴；

--hete：計(jì)算個(gè)體雜合率；

--out：指定輸出的文件路徑和文件名膏蚓。

程序運(yùn)行效果如下：

基因組文件中800個(gè)個(gè)體在1000個(gè)基因座中a1a2等位基因雜合的頻率保存在demo.hete文件中瓢谢。

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

第一列為基因型個(gè)體ID驮瞧；第二列為該個(gè)體a1a2等位基因雜合的頻率氓扛。

HIBLUP還提供了進(jìn)行個(gè)體篩選的參數(shù)--keep和--remove。所需要輸入的文件格式如id.filter.txt所示论笔，為一行一個(gè)的個(gè)體名稱采郎。

--keep：輸入包含分析中需要使用的個(gè)體ID的文件，分析時(shí)將只對(duì)文件中包含的個(gè)體進(jìn)行計(jì)算狂魔；

--remove：輸入包含分析中需要?jiǎng)h除的個(gè)體ID的文件蒜埋，分析時(shí)將刪除文件中包含的個(gè)體，對(duì)其他個(gè)體進(jìn)行計(jì)算最楷。

3整份、修改線程數(shù)

為了達(dá)到更高的計(jì)算效率，提高資源利用率籽孙，HIBLUP會(huì)默認(rèn)獲取OpenMP環(huán)境變量中最大線程數(shù)用于分析烈评。由于沒有指定線程數(shù)，這里默認(rèn)使用32個(gè)線程進(jìn)行計(jì)算犯建。

如果需要指定線程數(shù)讲冠，只需要在計(jì)算時(shí)使用--threads參數(shù)指定用于分析的線程數(shù)。比如將上述計(jì)算基因型頻率的命令行代碼更換為：

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

就可以只使用10個(gè)線程進(jìn)行計(jì)算胎挎，但計(jì)算效率也可能會(huì)受到影響沟启。