GWAS哈溫平衡檢驗(yàn)的方法

「什么是哈溫平衡？」

?  

哈迪-溫伯格(Hardy-Weinberg)法則 哈迪-溫伯格（Hardy-Weinberg）法則是群體遺傳中最重要的原理，它解釋了繁殖如何影響群體的基因和基因型頻率。這個(gè)法則是用Hardy,G.H (英國數(shù)學(xué)家) 和Weinberg，W.（德國醫(yī)生）兩位學(xué)者的姓來命名的，他們于同一年（1908年）各自發(fā)現(xiàn)了這一法則。他們提出在一個(gè)不發(fā)生突變、遷移和選擇的無限大的隨機(jī)交配的群體中，基因頻率和基因型頻率將逐代保持不變。---百度百科

?  

「怎么做哈溫平衡檢驗(yàn)？」

?  

「卡方適合性檢驗(yàn)！」，一個(gè)群體是否符合這種狀況，即達(dá)到了遺傳平衡，也就是一對等位基因的3種基因型的比例分布符合公式：p2+2pq+q2=1,p+q=1,(p+q)2=1.基因型MM的頻率為p2,NN的頻率為q2,MN的頻率為2pq。MN:MN：NN＝P2：2pq：q2。MN這對基因在群體中達(dá)此狀態(tài)，就是達(dá)到了遺傳平衡。如果沒有達(dá)到這個(gè)狀態(tài)，就是一個(gè)遺傳不平衡的群體。但隨著群體中的隨機(jī)交配，將會保持這個(gè)基因頻率和基因型分布比例，而較易達(dá)到遺傳平衡狀態(tài)。應(yīng)用Hardy-Weinberg遺傳平衡吻合度檢驗(yàn)方法，把計(jì)算得到的基因頻率代入，計(jì)算基因型平衡頻率，再乘以總?cè)藬?shù)，求得預(yù)期值（e）。把觀察數(shù)（O）與預(yù)期值（e）作比較，進(jìn)行χ2檢驗(yàn)。病例組和對照組的基因型分布的觀察值和預(yù)期值差異無顯著性（P>0.05），符合遺傳平衡定律.

?  

「哈溫平衡過濾和MAF過濾的區(qū)別？」

?  

之前，我對這兩個(gè)概念有點(diǎn)混淆，后來明白過來了。這兩個(gè)概念一個(gè)是對基因頻率進(jìn)行的篩選，一個(gè)是對基因型頻率進(jìn)行的篩選。對于一個(gè)位點(diǎn)“AA AT TT”，其中A的頻率為基因頻率，AA為基因型頻率。MAF直接是對基因頻率進(jìn)行篩選，而哈溫平衡檢驗(yàn)，則是根據(jù)基因型推斷出理想的（AA，AT，TT）的分布，然后和實(shí)際觀察的進(jìn)行適合性檢驗(yàn)，然后得到P值，根據(jù)P值進(jìn)行篩選。即P值越小，說明該位點(diǎn)越不符合哈溫平衡。

?  

「兩個(gè)目的：」

• 計(jì)算所有位點(diǎn)的哈溫檢測結(jié)果
• 刪除SNP中不符合哈溫平衡的位點(diǎn) 

1. 計(jì)算所有位點(diǎn)的HWE的P值

plink --bfile HapMap_3_r3_8 --hardy 

plink.hwe的數(shù)據(jù)格式：

• CHR 染色體
• SNP SNP的ID
• TEST 類型
• A1 minor 位點(diǎn)
• A2 major 位點(diǎn)
• GENO 基因型分布：A1A1, A1A2, A2A2
• O(HET) 觀測雜合度頻率
• E(HET) 期望雜合度頻率
• P 哈溫平衡的卡方檢驗(yàn)P-value值

結(jié)果預(yù)覽： 

2. 提取哈溫p值小于0.0001的位點(diǎn)

這里我們使用awk：

awk '{if($9 < 0.0001) print $0}' plink.hwe >plinkzoomhwe.hwe 

共有123個(gè)位點(diǎn)，其中UNAFF為45個(gè)位點(diǎn)。 

3.  設(shè)定過濾標(biāo)準(zhǔn)1e-4

plink --bfile HapMap_3_r3_8 --hwe 1e-4 --make-bed --out HapMap_3_r3_9
 

日志：

Options in effect:
 --bfile HapMap_3_r3_8
 --hwe 1e-4
 --make-bed
 --out HapMap_3_r3_9

515185 MB RAM detected; reserving 257592 MB for main workspace.
1073788 variants loaded from .bim file.
163 people (79 males, 84 females) loaded from .fam.
112 phenotype values loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 112 founders and 51 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.998136.
--hwe: 45 variants removed due to Hardy-Weinberg exact test.
1073743 variants and 163 people pass filters and QC.
Among remaining phenotypes, 56 are cases and 56 are controls.  (51 phenotypes
are missing.)
--make-bed to HapMap_3_r3_9.bed + HapMap_3_r3_9.bim + HapMap_3_r3_9.fam ...
done.
 

可以看到，共有45個(gè)SNP根據(jù)哈溫的P值過濾掉了，和上面手動計(jì)算的一樣。

4. 可視化

R代碼：

hwe<-read.table (file="plink.hwe", header=TRUE)
pdf("histhwe.pdf")
hist(hwe[,9],main="Histogram HWE")
dev.off()

hwe_zoom<-read.table (file="plinkzoomhwe.hwe", header=TRUE)
pdf("histhwe_below_theshold.pdf")
hist(hwe_zoom[,9],main="Histogram HWE: strongly deviating SNPs only")
dev.off()
 

哈溫的P值直方圖：

過濾掉SNP位點(diǎn)的P值： 

過濾后的結(jié)果文件

HapMap_3_r3_9.bed  HapMap_3_r3_9.bim  HapMap_3_r3_9.fam  HapMap_3_r3_9.log