【GWAS】使用PLINK評估族群分層的現象 @ 凡走過必留下痕跡~~~

在全基因體關聯性研究(Genome-Wide Association Study, GWAS)的分析中
干擾因子(confounder)的控制是重要的課題，其中最重要的干擾因子為族群分層(population stratification)
性狀(trait)或疾病的分布在不同族群中會不一樣，SNP的頻率在不同族群的分布也會不一樣
因此在探討兩者的關係時，若沒有考慮到族群分層，可能會出現偽陽性或偽陰性的結果
即使分析的研究對象都是來自單一族群，仍可能存在次分群的狀況
所以在完成QC後，進行相關性檢定前，需對研究族群進行主成分分析(Principle Component Analysis, PCA)
隨後將PCA結果作為共變數(covariate)，加入迴歸模型中以控制族群分層的影響

如同判斷親緣關係的流程，先篩選出獨立的SNPs，再進行PCA的計算
因此第一個步驟的PLINK指令為：

plink --bfile qc-final --exclude range high-LD-regions.txt --indep-pairwise 200 5 0.2 --out indep

--indep-pairwise代表用LD的方式來篩選SNPs，後面接的三個數字是定義篩選的條件
意思是取200個SNPs為一個視窗(window size)，計算彼此的r2，並以0.2為切點判斷是否有LD
此視窗計算完成後移動5個SNPs，再取200個SNPs為一個視窗重複一樣的計算，直到整條染色體算完
在基因組中，某些區域具有long range LD的特性，例如HLA，需額外排除才能篩選出真正互相獨立的SNPs
此步驟會產出兩個檔案：indep.prune.in是獨立的SNPs清單；indep.prune.out是具有LD的SNPs清單

接著利用互為獨立的SNPs來進行PCA

plink --bfile qc-final --extract indep.prune.in --pca header --out pca

PLINK的PCA分析，預設是產出20個PC值，因此pca.eigenvec檔內是每一位個案前20個PC的數值
接著用R繪製前10個PC值倆倆的散佈圖，並依個案所來自的族群標示不同顏色

pc <- read.table("pca.eigenvec", header = TRUE)
pop <- read.table("population.txt", header = TRUE)
pc <- merge(pc, pop, by = c("FID", "IID"))
color <- c("red", "green", "blue")
pairs(pc[3:12], col = color[pc$Population])

此次分析的資料是來自HapMap 3的亞洲人數據
因此在PC1對PC2的散佈圖中可以看出來研究對象被分成三個clusters
其中以來自CHB, CHD的cluster(紅點，綠點)彼此距離較近，來自JPT的cluster(藍點)距離另兩者都較遠
反應出這組資料確實是有群族分層的現象存在，在迴歸分析時應將這些PCs加入模式中
至於要加入幾個PCs？則視族群分層的狀況而定，並沒有一個定值

透過PCA方法，能夠幫忙檢視研究對象裡族群分層的狀況，也能協助判斷是否有離群值(outlier)存在
這個需將資料先與公開的資料庫合併，如：HapMap, 1000 Genomes等，再進行PCA
然後在PCA圖中觀察是否有個案遠離大部分的研究對象，或是偏離目標族群
利用這種定位方法，可以讓研究者獲得更同質的資料，也算是一種QC方法，如下圖的例子