「群体遗传学实战」第二课: 画出和文章几乎一样的PCA图

2020-04-15 本文已影响0人 xuzhougeng

主成分分析(PCA)是一种线性降维方法，能从纷繁复杂的数据中抽离出关键因素，用来区分不同的样本。这里我们不谈PCA背后的数学原理，只谈哪些软件能够处理数据，我找到了以下三款

无论使用哪款软件，始终都要记得它们最初为人类基因组设计，因此一定要关注和染色体有关的参数。例如Plink一定要加上 --allow-extra-chr允许非标准的染色体编号

首先，我们需要将VCF转成Plink格式. 由于VCF只包括样本编号，没有记录样本对应的群体信息，我们需要在转换的时候加入该信息。如果你的样本编号以来源信息+样本信息的方法是命名，那你可以用--id-delim，考虑到我们可能没有那么良好的命名习惯，所以这里直接用 --const-fid,

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

这一步会生成以.map,.nosex和.ped结尾的三个文件，这里是 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结尾的两个文件，这里就是watermelon_414acc.bim和watermelon_414acc.bed.

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

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

这一步会得到两个文件，一个是以.eigenval结尾的文件，记录特征值，用来计算每个PC所占的比重。另一个是.eigenvec结尾的文件，记录特征向量，用于坐标轴。

尽管你也可以使用GCTA计算PCA，但是他要求染色体的编号一定得是数字，不然在读取bim结尾的文件时，一定会报illegal chr number错误。此外，Plink在文档中提到, 如果想要自动离群值移除、大数据集处理，LD拟合，可以尝试用EIGENSOFT。

最终的结果，我们使用ggplot2进行可视化展示，我们得到了和原文一摸一样的图

8品种PCA

不难发现CN品种(只有一个样本，来自于Namibia),相对于其他品种都特别的远，就导致其他品种都挤到了一起，因此正文的PCA就剔除了CN这个品种。同时为了展示出C. lanatus的差异性，还单独分析C. mucosospermus和C. lanatus.

那么问题来了，我们如何在用Plink分析PCA时只选择部分样本呢? 我们可以使用Plink的--keep和--remove参数，只分析给定的样本。它们要求输入文件为两列，一列是样本所在的群体编号，一列是样本编号。以剔除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图会和系统发育树以及群体结构一起解释，相互验证。

此处是R代码，用到的数据可以从腾讯微云下载，链接：https://share.weiyun.com/5eFjaTZ 密码：3q514l