群體結構分析三種常用方法 (上篇)

寫在前面

在群體遺傳學和進化生物學相關的項目中，群體結構分析是最常見也是最初步的分析內容，可以幫助我們確認樣本分群是否符合預期以及檢測離群樣本。群體結構分析最常用的三種方法就是PCA、系統發生樹和祖先成分堆疊圖，下面我們將使用發表在Genome Rearch上的Gou et al,2014中的數據（60隻狗全基因組SNP）逐一講解,分為上下兩篇。

一、 PCA分析1. 簡單介紹PCA原理2. PCA實踐GCTA進行PCA分析EIGENSOFT進行PCA分析3. 繪製PCA散點圖小結參考文獻

一、 PCA分析1. 簡單介紹PCA原理

PCA (Principal Component Analysis) ，即主成分分析方法，是一種使用最廣泛的數據降維算法，通過正交變換將一組數量龐大且可能存在相關性的變量轉換為一組低維的線性不相關的變量。
對於一個全基因組測序樣本call出的SNP數目通常是百萬級甚至是千萬級別，比如本文所用數據為17M SNP，如果直接使用這些SNP位點作為指標對個體進行區分，就會由於信息過於冗雜而無法把握重點，並造成計算資源和時間的浪費。PCA分析的過程就是從千萬級別的SNP位點中提取關鍵信息，以便我們使用更少的變量就可以對樣本進行有效的刻畫和區分。

2. PCA實踐

在Gou et al, 2014文章中，作者使用GCTA進行PCA分析，在本文中，還會講解另外一種主流的PCA分析軟體--EIGENSOFT中的smartpca。

軟體下載連結：
GCTA
EIGENSOFT

GCTA進行PCA分析

首先準備GCTA中PCA模塊所需的輸入文件，GCTA可以直接讀取.bed , .bim , .fam文件，使用plink將vcf格式文件轉換成上述三種文件,同時進行簡單的過濾(--geno 0.05)。

#!/bin/bash

plink=/software/biosoft/software/plink/plink
time $plink --vcf 60_dog.merge.vcf --geno 0.05 --dog --make-bed --out 60_dog_geno0.05 && echo "---- vcf2bed done "

然後使用GCTA軟體，--make-grm 生成個體對之間的遺傳關係矩陣(genetic relationship matrix),並將GRM的下三角元素保存為二進位文件.grm.id , .grm.bin ， .grm.N.bin。
這裡要強調一點，對於非人物種，要設定好染色體數目參數，--autosome-num。否則會因為不識別26以上的染色體編號而報錯。

gcta=/software/biosoft/software/gcta_1.92.1beta6/gcta64
time $gcta --bfile 60_dog_geno0.05 --make-grm --autosome-num 38 --out 60_dog_geno0.05 && echo "----make-grm done----"

最後就是要進行PCA分析
使用 --pca 設置要生成主成分的數目，這裡設置為4，一般來說就可以刻畫出群體結構。
這一步會生成 .eigenval 和 .eigenvec 兩個文件。.eigenval文件為各主成分可解釋遺傳信息的比例，.eigenvec文件為每個樣本在top4主成分上的分解值。

time $gcta --grm 60_dog_geno0.05 --pca 4 --out 60_dog_geno0.05 && echo "----gcta-pca done----"

補充樣本所屬子群信息
在生成的.eigenvec文件中，並沒有包含樣本所屬子群信息，在進行可視化之前，需要添加上樣本群體信息。樣本較少時可以手動添加，樣本數目較多時推薦寫代碼完成，避免錯誤。文章Gou et al,2014中根據採樣地點劃分子群。添加上子群信息後，.eigenvec長這個樣子。

EIGENSOFT進行PCA分析

使用plink將vcf格式文件轉換成.bed , .bim , .fam文件，步驟同上，這裡不再贅述。

smartpca需要提供三種輸入文件，文件1記錄每個樣本在每個SNP位點的信息，文件2記錄每個SNP位點的信息，文件3記錄每個樣本的性別和所屬子群信息。
熟悉plink格式的同學應該知道.bed文件儲存的就是genotype信息，.bim文件存儲的就是SNP位點信息，那麼我們只需要自己生成文件3就可以了。文件3格式如下，第一列為樣本ID，第二列為性別（M，F，U），第三列為所屬子群。

接下來把所有參數寫入一個文件中，就可以運行smartpca了，其中 numoutevec 是輸出主成分的數目。

運行smartpca，生成 .evec 和 .eval 文件。

time /software/biosoft/software/EIG-6.1.4/bin/smartpca -p PCA.par > smartpca.log && echo "----smartpca----done"

截圖來看下 .evec 文件，兩種軟體運行結果差不多。

3. 繪製PCA散點圖

通過得到的.evec和eval兩個文件，這裡使用python完成可視化，代碼如下：

import matplotlib.pyplot as plt
import collections
import re

hash = collections.OrderedDict()
eval_file = open("60_dog_geno0.05.eigenval","r")
evec_file = open("60_dog_geno0.05.eigenvec","r")

###  從.eval文件中讀取top4主成分
eval_1 = eval_file.readline()
eval_2 = eval_file.readline()
eval_3 = eval_file.readline()
eval_4 = eval_file.readline()


###  從.evec文件中讀取在pc上的值
for i in evec_file:
    if re.match(r'#',i.strip()):
        continue
    tmp = i.strip().split()[1:]
    if tmp[-1] in hash.keys():
        hash[tmp[-1]]['1st'].append(eval(tmp[1]))
        hash[tmp[-1]]['2nd'].append(eval(tmp[2]))
        hash[tmp[-1]]['3rd'].append(eval(tmp[3]))
        hash[tmp[-1]]['4th'].append(eval(tmp[4]))
    else:
        hash[tmp[-1]] = hash.get(tmp[-1],collections.OrderedDict())
        hash[tmp[-1]]['1st'] = hash[tmp[-1]].get('1st',[])
        hash[tmp[-1]]['2nd'] = hash[tmp[-1]].get('2nd',[])
        hash[tmp[-1]]['3rd'] = hash[tmp[-1]].get('3rd',[])
        hash[tmp[-1]]['4th'] = hash[tmp[-1]].get('4th',[])
        hash[tmp[-1]]['1st'].append(eval(tmp[1]))
        hash[tmp[-1]]['2nd'].append(eval(tmp[2]))
        hash[tmp[-1]]['3rd'].append(eval(tmp[3]))
        hash[tmp[-1]]['4th'].append(eval(tmp[4]))

### 繪製散點圖，這裡只繪製pc1和pc2,其他pc可按照下面代碼逐一畫出。
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(1, 1, 1)            ### 2D figure
#ax = fig.add_subplot(111,projection='3d')              ### 3D figure
mark = ['o','v','s','*','x','+']*10         ### 設置散點性狀
col = ['b','g','y','r','c']*10               ### 設置散點顏色
for m,n in enumerate(hash.keys()):    ### 逐點畫出
    ax.scatter(hash[n]['1st'],hash[n]['2nd'],c=col[m],s=75,marker=mark[m],label=n,alpha=0.7)
ax.legend(loc='best',scatterpoints=1,fontsize='12',framealpha=0)
ax.set_xlabel('Eigenvector1 ({}%)'.format(float('%.2f' % eval(eval_1))),fontsize=13,fontweight='bold')
ax.set_ylabel('Eigenvector2 ({}%)'.format(float('%.2f' % eval(eval_2))),fontsize=13,fontweight='bold')

#### 保存圖片
save_FileName = '60dog_geno0.05_gcta_pc1_pc2.png'
plt.savefig(save_FileName,dpi=400,bbox_inches='tight')

PCA散點圖如下

小結

PCA運用方差分解，將海量SNP的差異反映在二維坐標圖上，坐標軸軸取對樣品差異性解釋度最高的兩個或三個特徵值，樣本SNP位點信息越相似，反映在PCA圖中的距離越近。PCA可以幫助我們清楚掌握手上樣本的群體結構，並有效檢測出離群樣本，為下遊分析減少不必要的麻煩，同時，實現PCA分析的方法很多，PCA結果圖也簡單易懂，稱得上是低調實力派。
最後強調一點，PCA分析是降維，不是聚類，希望大家不要搞混了。
對PCA理論和分析過程感興趣的同學，可以進一步看下這篇文章。

謝謝大家賞臉看到這裡，在本篇文章中，主要學習了PCA分析的基本過程，在下篇文章中，我們將接著學習系統發生樹構建和祖先成分分析，精彩不容錯過，請持續關注我們！

參考文獻

Genome Res. 2014 Aug;24(8):1308-15. doi: 10.1101/gr.171876.113. Epub 2014 Apr 10.

http://blog.csdn.net/zhongkelee/article/details/44064401