主成分分析(PCA)是一種線性降維方法，能從紛繁復(fù)雜的數(shù)據(jù)中抽離出關(guān)鍵因素，用來區(qū)分不同的樣本。這里我們不談PCA背后的數(shù)學(xué)原理，只談哪些軟件能夠處理數(shù)據(jù)，我找到了以下三款

• Plink: https://www.cog-genomics.org/plink/
• GCTA: https://cnsgenomics.com/software/gcta/
• EIGENSOFT: https://github.com/DReichLab/EIG

無論使用哪款軟件，始終都要記得它們最初為人類基因組設(shè)計，因此一定要關(guān)注和染色體有關(guān)的參數(shù)。例如Plink一定要加上 --allow-extra-chr允許非標(biāo)準(zhǔn)的染色體編號

首先，我們需要將VCF轉(zhuǎn)成Plink格式. 由于VCF只包括樣本編號，沒有記錄樣本對應(yīng)的群體信息，我們需要在轉(zhuǎn)換的時候加入該信息。如果你的樣本編號以來源信息+樣本信息的方法是命名，那你可以用--id-delim，考慮到我們可能沒有那么良好的命名習(xí)慣，所以這里直接用 --const-fid,

plink --vcf watermelon_414acc_SNP2.vcf.gz --recode --out watermelon_414acc --const-fid --allow-extra-chr

這一步會生成以.map,.nosex和.ped結(jié)尾的三個文件，這里是 watermelon_414acc.map, watermelon_414acc.nosex, watermelon_414acc.ped

之后我們還需基于PED生成一個bed文件，這里bed指的是binary ped，并非基因組瀏覽器用到BED文件。

plink --allow-extra-chr --file watermelon_414acc --noweb --make-bed --out watermelon_414acc

這一步會生成以.bim,bed結(jié)尾的兩個文件，這里就是watermelon_414acc.bim和watermelon_414acc.bed.

最后就可以計算利用Plink計算PCA

plink --allow-extra-chr --threads 20 -bfile watermelon_414acc --pca 20 --out watermelon_414acc

這一步會得到兩個文件，一個是以.eigenval結(jié)尾的文件，記錄特征值，用來計算每個PC所占的比重。另一個是.eigenvec結(jié)尾的文件，記錄特征向量，用于坐標(biāo)軸。

盡管你也可以使用GCTA計算PCA，但是他要求染色體的編號一定得是數(shù)字，不然在讀取bim結(jié)尾的文件時，一定會報illegal chr number錯誤。此外，Plink在文檔中提到, 如果想要自動離群值移除、大數(shù)據(jù)集處理，LD擬合，可以嘗試用EIGENSOFT。

最終的結(jié)果，我們使用ggplot2進行可視化展示，我們得到了和原文一摸一樣的圖

8品種PCA

不難發(fā)現(xiàn)CN品種(只有一個樣本，來自于Namibia),相對于其他品種都特別的遠(yuǎn)，就導(dǎo)致其他品種都擠到了一起，因此正文的PCA就剔除了CN這個品種。同時為了展示出C. lanatus的差異性，還單獨分析C. mucosospermus和C. lanatus.

那么問題來了，我們?nèi)绾卧谟肞link分析PCA時只選擇部分樣本呢? 我們可以使用Plink的--keep和--remove參數(shù)，只分析給定的樣本。它們要求輸入文件為兩列，一列是樣本所在的群體編號，一列是樣本編號。以剔除CN樣本為例，它的樣本編號為WM439，所在群體是0。

echo '0\tWM439' > remove.txt
plink --remove remove.txt --allow-extra-chr -bfile watermelon_414acc --pca 20 --out watermelon_414acc_no_cn

然后用相同的代碼畫圖，和文章的正文完全一致。

不包括CN品種

通常而言，PCA圖會和系統(tǒng)發(fā)育樹以及群體結(jié)構(gòu)一起解釋，相互驗證。

此處是R代碼，用到的數(shù)據(jù)可以從騰訊微云下載，鏈接：https://share.weiyun.com/5eFjaTZ 密碼：3q514l