PCA方法校正群體結構，GWAS該用多少個主成分？
該選擇多少個主成分
群體結構（population structure），或者說群體分層（population stratification），是由於個體之間非隨機交配而導致的群體中亞群之間等位基因頻率的系統差異。這種系統差異，是全基因組關聯研究（GWAS）中影響非常大的混淆變數，可以造成非常大的假陽性。
舉個簡單的模擬例子 [1]，當 GWAS 中不存在群體分層時，得到的結果會是比較真實可靠的：
當樣本存在一定程度的群體分層現象時，會出現一些假陽性訊號：
當群體分層現象非常嚴重時，bonferroni correction 校正也沒什麼用，大量位點都會超過 bonferroni correction 的閾值：
為了儘量降低群體結構的影響，通常會先對基因組進行主成分分析（PCA），然後在做 GWAS 時會加入主成分（principal components, PCs）作為協變數。
但問題就來了，該選擇多少個主成分去校正群體結構？PCA 個數的選擇對結果影響很大。如果選擇的個數太少，無法有效校正群體結構，假陽性仍然會很大。但如果選擇的個數太多，會影響 GWAS 的 power。 下面就說說常見的幾種方法。
直接選取前 k 個主成分
最簡單直接的方法就是人為選擇前 k 個 PCs 作為協變數，比如直接選取前 5 個或者前 10個。早期的文獻通常推薦使用前 10 個 PCs作為協變數，校正群體結構 [1]。
不過，這種方法過於簡單粗暴。在人群數量和樣本數量快速增長、一個 GWAS 能達到幾萬人甚至幾十萬人的今天， 這樣的粗暴方法往往並不足以校正群體結果。
所以，這種方法雖然簡單，但並不推薦。
基於 PCA 散點圖或者 ANOVA
如果要更為可靠地選取 PCs 數量，可以繪製用 eigenvector 繪製散點圖，選擇可以將群體有效分開前 k 個 的主成分。比如下面這張圖，前兩個 PCs 可以將 3 個群體分開，而 PC3、PC4 無法將三個群體分開。這時，選擇 PC1 和 PC2 作為 GWAS 的協變數就足夠了。如果強行把 PC3 和 PC4 也加進去，反而會帶來很大的 bias。
進行 PCA 並且畫圖的操作過程可看看
覺得畫圖不夠客觀，或者太過麻煩，如果知道群體的個數的話也可以對主成分進行 ANOVA 分析，檢驗主成分在不同人群中是否顯著，選擇顯著的前 k 個主成分。
twstats 方法（推薦）
第二種畫圖的方法觀察起來還是有些主觀，如果不能很好定義人群，ANOVA 的方法也不太好用。是否有更好的方法？
這裡推薦 Patterson 在 2006 年發表的 EIGENSTRAT 文章中的 twstats 方法（compute number of statistically significant principal components） [2]。它基於 Tracy–Widom statistics，對各個主成分進行顯著性檢驗。在模擬結果中，Tracy–Widom statistics 的顯著性檢驗結果與 ANOVA 比較吻合，可靠性不錯。
這種方法整合在 EIGENSOFT 的 twtable 中。Plink 或者 EIGENSTRAT 的 PCA 結果可直接用來計算：
# 利用 Plink 計算 PCA，輸出前 50 個主成分
plink --bfile yourfile --pca 50 --out yourfile_pca
# 下載解壓，如果下面的網址無法下載，可以試試 https://storage.googleapis.com/broad-alkesgroup-public/EIGENSOFT/EIG-6.1.4.tar.gz
wget https://data.broadinstitute.org/alkesgroup/EIGENSOFT/EIG-6.1.4.tar.gz
tar xzvf EIG-6.1.4.tar.gz
# 利用 twstats 進行顯著性檢驗
/bin/twstats -t twtable -i yourfile_pca.eigenval -o yourfile_pca_number
在結果輸出檔案中選擇 P < 0.05 的前 k 個主成分：
不過，用 twstats 評估顯著性時要注意：The twstats program assumes a random set of markers, and should not be used on data sets of ancestry-informative markers, as admixture-LD may violate its underlying assumptions.。也就是說，如果你的基因型資料都是 ancestry-informative markers，就不能用這種方法。
基於可解釋方差（推薦）
除了基於 Tracy–Widom statistics 檢驗主成分的顯著性外，還可以根據每個主成分的可解釋方差計算。一般，選取累計解釋 80-90% 的前 k 個主成分就足夠。Plink、GCTA 等工具不能輸出各個主成分的可解釋方差，要這個資訊的話可以用 vcfR、SNPRelate、bigsnpr、pcadapt 等 R 包。
SNPRelate 的平行計算速度比較快，以它為例，計算 PCA 並且得到可解釋方差：
# from shiyanhe and zhaozhuji.net
# 從 Bioconductor 安裝 SNPRelate 包和它依賴的 gdsfmt 包
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("gdsfmt")
BiocManager::install("SNPRelate")
# 載入 gdsfmt 和 SNPRelate 包
library(gdsfmt)
library(SNPRelate)
# 輸入 PLINK 檔案路徑
bed.fn <- "/your_folder/your_plink_file.bed"
fam.fn <- "/your_folder/your_plink_file.fam"
bim.fn <- "/your_folder/your_plink_file.bim"
# 將 PLINK 檔案轉為 GDS 檔案
snpgdsBED2GDS(bed.fn, fam.fn, bim.fn, "test.gds")
# 讀取 GDS 檔案
genofile <- snpgdsOpen("test.gds")
# 根據 LD 過濾 SNPs，閾值根據需要設定
set.seed(1000)
snpset <- snpgdsLDpruning(genofile, ld.threshold=0.2)
# 選擇 SNP pruning 後要保留的 SNP
snpset.id <- unlist(unname(snpset))
# 計算 PCA，num.thread 是並行的執行緒數
pca <- snpgdsPCA(genofile, snp.id=snpset.id, num.thread=10)
# 以百分比形式輸出 variance proportion
print(pca$varprop*100)
當然，也可以繪製碎石圖（ scree plot）來觀察：
# 繪製前 30 個主成分的碎石圖
# from shiyanhe and zhaozhuji.net
library(ggplot2)
K= 30
qplot(x = 1:K, y = (pca$varprop[1:K]), col = "red", xlab = "PC", ylab = "Proportion of explained variance") +
geom_line() + guides(colour = FALSE) +
ggtitle(paste("Scree Plot - K =", K))
這樣就可以畫出碎石圖：
總結
總的來說，這裡推薦第三種方法 twstats 和第四種基於累計可解釋方差的方法。在實際應用中，建議同時結合這兩種方法。首先用 twstats 方法計算各個主成分顯著性，再計算各個主成分的可解釋方差，然後選取 P 值顯著且累計可解釋方差在 80-90% 的前 k 個主成分。這樣可以讓結果更為可靠，選擇更恰當的主成分個數。
參考：
[1] Zhao H, Mitra N, Kanetsky P A, et al. A practical approach to adjusting for population stratification in genome-wide association studies: principal components and propensity scores (PCAPS)[J]. Statistical applications in genetics and molecular biology, 2018, 17(6).
[2] Patterson N, Price A L, Reich D. Population structure and eigenanalysis[J]. PLoS genet, 2006, 2(12): e190.
[3] https://www.biostars.org/p/126274/
[4] https://www.biostars.org/p/305469/
[5] https://www.biostars.org/p/335605/
[6] https://www.cnblogs.com/chenwenyan/
[7] twstats文件：https://github.com/chrchang/eigensoft/blob/master/POPGEN/README