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

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

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

无论使用哪款软件，始终都要记得它们最初为人类基因组设计，因此一定要关注和染色体有关的参数。例如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

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

「群体遗传学实战」第三课: 如何对SNP位点进行过滤
往期教程「群体遗传学实战」第一课: 对SNP位点进行注释「群体遗传学实战」第二课: 画出和文章几乎一样的PCA图 SNP过滤有两种情况,一种是仅根据位点质量信息(测序深度,回帖质量等)对SNP进行 ...
「群体遗传学实战」第一课: 对SNP位点进行注释
文章我们用于实战的数据集来自于2019年发表于NG的西瓜文章,它提供了GATK过滤后的SNP数据集用于分析,并且附录提供了完整的样本信息.该SNP数据集包括414个西瓜,2000万个SNP信息,数据 ...
计算机画大熊猫教学设计,第二课画大熊猫教学设计.docx
今泰学院推荐文档第 PAGE 第 PAGE 2 页共 NUMPAGES 2 页第二课画大熊猫教学设计Lesson 2 teaching design of drawing giant panda ...
pyhton 画出音频文件的波形图和频谱图
pyhton 画出音频文件的波形图和频谱图 # -*- coding:utf-8 -*- import wave import struct import numpy as np import mat ...
python画人民大会堂_太震撼了，我用python画出全北京的公交线路动图
原标题:太震撼了,我用python画出全北京的公交线路动图今天教大家用pyecharts制作北京市公交线路动态图,这应该是全网唯一一篇能正常运行的教程一.获取百度秘钥首先,本项目需要引用百度地图 ...
wps 甘特图_「WPS办公助手」3分钟画出高逼格流程图，工作学习规划一目了然，老板都要为你点赞...
今天小金想和大家唠唠 WPS 两个独特的办公工具:「WPS流程图」和「WPS脑图」.可能很多人潜意识里觉得流程图和脑图比较麻烦很难画,甚至还有一些盆友从来没有听过.但实际上无论工作学习还是日常生活,流 ...
「Vue实战」武装你的前端项目
1. 接口模块处理 1.1 axios二次封装很基础的部分,已封装好的请跳过.这里的封装是依据 JWT import axios from 'axios'import router from '.. ...
思维导图怎么画？怎么画出简单又漂亮的思维导图
思维导图怎么画?现在主要的两种方式:一种是手绘,一种是借助思维导图网站进行绘制.思维导图是一种思维整理工具,对于上班族来说,使用思维导图可以帮助你提高工作效率.对于学生党来说使用思维导图可以锻炼逻辑能 ...
怎么查询redis缓存的数据_阿里开发十年写出这份「Redis简明教程」+「Redis实战」请你查收...
Redis是啥?用Redis官方的话来说就是: Redis is an open source (BSD licensed), in-memory data structure store, used ...

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

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

最新文章

热门文章