GWAS全基因組關(guān)聯(lián)分析第二期：數(shù)據(jù)質(zhì)控

上一期我們了分享了GWAS分析需要的數(shù)據(jù)格式，以及不同格式之間的轉(zhuǎn)換。現(xiàn)在我們已經(jīng)準(zhǔn)備好了表型數(shù)據(jù)和基因數(shù)據(jù)，是不是就想馬上進(jìn)行關(guān)聯(lián)分析了？心急吃不了熱豆腐，為了提高關(guān)聯(lián)分析結(jié)果的準(zhǔn)確性，需要對(duì)數(shù)據(jù)進(jìn)行質(zhì)控，去掉不合格的樣本和變異數(shù)據(jù)。

1 SNP及個(gè)體缺失過(guò)濾 

人工采集的數(shù)據(jù)，可能存在位點(diǎn)基因型和個(gè)體基因數(shù)據(jù)缺失（表型缺失的直接去掉），這些缺失數(shù)據(jù)影響關(guān)聯(lián)分析的準(zhǔn)確性，需要將缺失率控制在一定標(biāo)準(zhǔn)以下。建議首先以寬松的閾值（0.2；> 20%）過(guò)濾SNP和個(gè)體，從而過(guò)濾掉缺失程度很高的SNP和個(gè)體；再使用更嚴(yán)格的閾值過(guò)濾（(0.02；> 2%）。

# SNP缺失過(guò)濾
$plink --noweb --bfile $project.raw.mark --geno 0.2 --allow-no-sex --make-bed --out ${project}.filter.mds1

# 個(gè)體缺失過(guò)濾
$plink --noweb --bfile ${project}.filter.mds1 --mind 0.2 --allow-no-sex --make-bed --out ${project}.filter.mds2

注意：以上步驟更換更嚴(yán)格的參數(shù)再過(guò)濾一遍。

2 性別和親緣關(guān)系檢測(cè)（可選） 

性別檢測(cè)基于X染色體近交系（純合子性）估計(jì)，一般女性受試者的F值 < 0.2，男性受試者的F值 > 0.8，不滿足這些要求的被標(biāo)記為“PROBLEM”。

# 性別檢測(cè)
$plink --noweb --bfile ${project}.raw.mark --check-sex

# 輸出結(jié)果保存在plink.sexcheck文件中，提取性別異常個(gè)體
$grep "PROBLEM" plink.sexcheck | awk '{print $1,$2}' >sex_removelist.txt

# 刪除性別異常個(gè)體（不建議刪除，除非明確該樣本數(shù)據(jù)有污染）
$plink --noweb --bfile ${project}.raw.mark --remove sex_removelist.txt --make-bed --out ${project}.raw.mark2

親緣關(guān)系檢測(cè)基于遺傳信息，判斷樣本親緣關(guān)系的指標(biāo)分為狀態(tài)同源（identical by state，IBS）和血緣同源（Identity By Descent，IBD），通常IBD無(wú)法直接觀察，但I(xiàn)BS可以通過(guò)兩個(gè)個(gè)體基因型算出（如下圖），再根據(jù)IBS以及等位基因頻率的分布推斷IBD。

 
# 親緣關(guān)系檢測(cè)
$plink --noweb --bfile ${project}.raw.mark --genome

# 輸出文件保存在plink.genome文件中，提取親緣關(guān)系異常的樣本
sed 's/^\s\+//' plink.genome | sed 's/\s\+/\t/g' | awk -v dst=0.85 'NR>2 {if($12 > dst) {print $1,$2; print $3,$4}}' | sort | uniq >genome_removelist.txt

# 刪除親緣關(guān)系異常個(gè)體（不建議刪除）
$plink --noweb --bfile ${project}.raw.mark --remove genome_removelist.txt --make-bed --out ${project}.raw.mark2

3 哈溫平衡過(guò)濾 

哈迪-溫伯格（Hardy-Weinberg）法則是群體遺傳中最重要的原理，提出在一個(gè)不發(fā)生突變、遷移和選擇的無(wú)限大的隨機(jī)交配的群體中（理想狀態(tài)下），基因頻率和基因型頻率將逐代保持不變。一對(duì)等位基因的3種基因型分布比例符合以下規(guī)律：
(p + q)^2 = 1 等價(jià)于 p^2 + 2pq + q^2 = 1
注：p和q分別表示兩個(gè)等位基因頻率，且p + q = 1。
$plink --noweb --bfile ${project}.raw.mark --hwe 1e-10 --hwe-all --make-bed --out ${project}.filter.haw

4 最小等位基因頻率過(guò)濾 

最小等位基因頻率（MAF）通常是指在給定人群中的不常見(jiàn)的等位基因發(fā)生頻率。

MAF如果非常小，比如低于0.02，那么意味著大部分位點(diǎn)都是相同的基因型，這些位點(diǎn)貢獻(xiàn)的信息非常少，增加假陽(yáng)性；更有甚者M(jìn)AF為0，即所有位點(diǎn)只有一種基因型，這些位點(diǎn)沒(méi)有貢獻(xiàn)信息，放在計(jì)算中增加計(jì)算量，沒(méi)有意義，所以要根據(jù)MAF進(jìn)行過(guò)濾。

# 最小等位基因頻率過(guò)濾（這里MAF閾值設(shè)為0.05）
$plink --noweb --bfile ${project}.raw.mark --maf 0.05 --allow-no-sex --make-bed --out ${project}.filter.maf

5 群體分層 

群體分層（Population stratification）：是最常見(jiàn)的差異來(lái)源，指的是case/control組的樣本來(lái)自于不同的祖先群體，其分型結(jié)果自然是有差異的。

不同群體SNP頻率不一樣，導(dǎo)致后面做關(guān)聯(lián)分析的時(shí)候可能出現(xiàn)假陽(yáng)性位點(diǎn)（不一定是顯著信號(hào)位點(diǎn)與該表型有關(guān)，可能是與群體SNP頻率差異有關(guān)），因此我們需要在關(guān)聯(lián)分析前對(duì)群體分層校正。

# 主成分分析
$plink --noweb --bfile ${project}.raw.mark --pca 10 --out pca

# 提取離群樣本
根據(jù)主成分分析結(jié)果，繪圖展示，確定離群樣本，寫(xiě)入pca_removelist.txt文件

# 刪除離群個(gè)體（可選）
$plink --noweb --bfile ${project}.raw.mark --remove pca_removelist.txt --make-bed --out ${project}.filter.pc

6 雜合性過(guò)濾 

雜合性是指某一個(gè)位點(diǎn)上含有一對(duì)及其以上的不同的等位基因。包括同系合性和同種合性。群體遺傳多態(tài)性的均勻度的度量常采用雜合度作為參數(shù)。雜合性是在同源染色體上的一個(gè)或多個(gè)位點(diǎn)上有不同等位基因存在的狀態(tài)，是種群的基本屬性之一。

# 連鎖過(guò)濾（LD），得到不連鎖的SNP
$plink --noweb --bfile ${project}.raw.mark --indep-pairwise 50 5 0.2 --out indepSNP

# 提取不連鎖的SNP進(jìn)行雜合性分析
$plink --noweb --bfile ${project}.raw.mark --extract indepSNP.prune.in --het --out hetSNP

# 提取雜合度較高的個(gè)體
sed 's/^\s\+//' hetSNP.het | sed 's/\s\+/\t/g' | awk -v f=0.35 'NR>1 {if(($5-$3)/$5 > f) {print $1,$2}}' >hetSNP_removelist.txt

# 刪除雜合度高的個(gè)體（可選）
$plink --noweb --bfile ${project}.raw.mark --remove hetSNP_removelist.txt --make-bed --out ${project}.filter.het

以上就是本期分享的內(nèi)容，下一期我們將講解GWAS關(guān)聯(lián)分析。