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

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

1 SNP及個體缺失過濾 

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

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

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

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

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

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

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

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

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

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

 
# 親緣關(guān)系檢測
$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)系異常個體（不建議刪除）
$plink --noweb --bfile ${project}.raw.mark --remove genome_removelist.txt --make-bed --out ${project}.raw.mark2

3 哈溫平衡過濾 

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

4 最小等位基因頻率過濾 

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

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

# 最小等位基因頻率過濾（這里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）：是最常見的差異來源，指的是case/control組的樣本來自于不同的祖先群體，其分型結(jié)果自然是有差異的。

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

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

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

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

6 雜合性過濾 

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

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

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

# 提取雜合度較高的個體
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

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

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