使用基因型数据计算近交系数

 

001、plink --het （纯合度近交系数）

root@PC1:/home/test# ls
outcome.map  outcome.ped
root@PC1:/home/test# plink --file outcome --het 1> /dev/null
root@PC1:/home/test# ls
outcome.map  outcome.ped  plink.het  plink.log
root@PC1:/home/test# head plink.het      ## 最后一列F为近交系数
    FID     IID       O(HOM)       E(HOM)        N(NM)            F
   DOR1    DOR1        28943    2.669e+04        42629       0.1412
   DOR2    DOR2        29737    2.669e+04        42629       0.1911
   DOR3    DOR3        28395    2.669e+04        42629       0.1068
   DOR4    DOR4        28430    2.669e+04        42629        0.109
   DOR5    DOR5        29235    2.669e+04        42629       0.1596
   DOR6    DOR6        28082    2.669e+04        42629       0.0872
   DOR7    DOR7        27779    2.669e+04        42629      0.06819
   DOR8    DOR8        29014    2.669e+04        42629       0.1457
   DOR9    DOR9        28135    2.669e+04        42629      0.09053

 

002、ROH 近交系数

root@PC1:/home/test# ls
outcome.map  outcome.ped                                                        ## plink检测ROH
root@PC1:/home/test# plink --file outcome --homozyg-snp 30 --homozyg-kb 1000 --homozyg-density 1000 --homozyg-gap 1000 --homozyg-window-snp 50 --homozyg-window-het 1 --homozyg-window-missing 1 --out roh 1> /dev/null
root@PC1:/home/test# ls
outcome.map  outcome.ped  roh.hom  roh.hom.indiv  roh.hom.summary  roh.log                                 ## 统计染色体的长度
root@PC1:/home/test# cut -f 1 outcome.map | sort -k 1n,1n| uniq | while read i; do grep -w "^$i" outcome.map | sort -k 4rn,4rn | head -n 1 >>chr_len.txt; done
root@PC1:/home/test# ls
chr_len.txt  outcome.map  outcome.ped  roh.hom  roh.hom.indiv  roh.hom.summary  roh.log
root@PC1:/home/test# geno_len=$(awk '{sum += $4} END {print sum / 1000}' chr_len.txt )  ## 计算基因组的长度                                           
root@PC1:/home/test# awk -v a=$geno_len '{print $2, $(NF - 1)/a}' roh.hom.indiv | head       ## 计算基因组近交系数
IID 0
DOR1 0.141861
DOR2 0.193722
DOR3 0.115116
DOR4 0.114285
DOR5 0.162598
DOR6 0.0992099
DOR7 0.0779776
DOR8 0.14785
DOR9 0.0968709

 

003、gcta计算近交系数

root@PC1:/home/test# ls
outcome.map  outcome.ped
root@PC1:/home/test# plink --file outcome --make-bed 1> /dev/null
root@PC1:/home/test# ls
outcome.map  outcome.ped  plink.bed  plink.bim  plink.fam  plink.log
root@PC1:/home/test# /home/software/gcta_v1.94.0Beta_linux_kernel_3_x86_64/gcta_v1.94.0Beta_linux_kernel_3_x86_64_static --bfile plink --ibc --out result 1> /dev/null
root@PC1:/home/test# ls
outcome.map  outcome.ped  plink.bed  plink.bim  plink.fam  plink.log  result.ibc  result.log
root@PC1:/home/test# head result.ibc
FID     IID     NOMISS  Fhat1   Fhat2   Fhat3
DOR1    DOR1    42629   0.0462481       0.154288        0.100268
DOR2    DOR2    42629   0.0309685       0.217734        0.124351
DOR3    DOR3    42629   0.0261114       0.114167        0.0701392
DOR4    DOR4    42629   -0.0207756      0.130706        0.0549651
DOR5    DOR5    42629   -0.0187196      0.189537        0.0854086
DOR6    DOR6    42629   -0.0154794      0.106942        0.0457311
DOR7    DOR7    42629   -0.0635135      0.0992687       0.0178776
DOR8    DOR8    42629   -0.0352431      0.177956        0.0713562
DOR9    DOR9    42629   0.0317887       0.0974181       0.0646034