群体结构 |
您所在的位置:网站首页 › pca图的意义 › 群体结构 |
概念
PCA(principal components analysis)即主成分分析。主成分分析也称主分量分析,旨在利用降维的思想,把多指标转化为少数几个综合指标。 在统计学中,主成分分析PCA是一种简化数据集的技术。它是一个线性变换的过程。这个变换把数据变换到一个新的坐标系统中,使得任何数据投影的第一大方差在第一个坐标(第一主成分)上,第二大方差在第二个坐标(第二主成分)上,依次类推。主成分分析经常用于减少数据集的维数,同时保持数据集的对方差贡献最大的特征。 PCA图形的解读PCA分析虽然朴实无华,其实就是散点图,但它应用十分广泛,能够帮我们解决很多生物学问题。PCA分析应用情境是:在某些情况下,生物数据实在过于复杂。 例如:对一个群体进行重测序,得到的SNP位点数可能是百万级别的。如果我们直接使用百万级别的SNP信息作为指标对个体进行区分,就会显得信息过于庞大而无法把握重点。PCA分析过程就是从这百万级别的信息中提取关键的信息,以便我们使用更少的标记就可以对样本进行有效区分。这些被提取出的信息,按照其效应从大到小排列,我们称之为主成分1(Principal component1)、主成分2、主成分3… … 运用1:群体结构分析在实际文章中,我们不仅仅使用PC1和PC2来对样本群体进行区分。从数学上理解,PCA分析的过程就是从大量数据指标中提取关键信息的过程。但PC1或PC2对总体信息的解释程度总是有限的。我们将之称为PCn对总体方差解释的百分比。一般重测序的PCA分析结果中,PC1对总体信息的解释比例在3~10%之间。所以,我们也需要关注一下其他主成分的分类效果。 例如在家蚕重测序文章中,分别使用主成分1和2绘图(左图)以及主成分3与主成分4绘图(右图)。两个聚类结果呈现了不同的意义。在PC1和PC2的聚类图中,将野生蚕和家蚕区分开了两个群体。而在PC3和PC4的聚类中,则分离出了两个来自江南地区高产丝量的品种。 家蚕PCA分析结果所以,从生物学层面理解,PCA分析的过程就是信息浓缩的过程,会从原始的各个SNP位点信息中提取相似的信息,浓缩为新的变量PC1、PC2、PC3…. 输出。所以不同的主成分可能会对应不同的生物学意义,产生不同的聚类分类效果。 运用2:检测离群样本例如,在上图(右)中,两个高产的品种就属于离群样本。如果你材料已知都是来源同一品种的个体,这种离群样本可能就意味着在采样或测序过程中,出现了样本混淆。如果这些材料后续用于GWAS分析,个别样本出现离群则考虑要把这些离群样本剔除。当然,如果大量样本离群或出现群体分层(例如,上图的左图,明显分层为两个亚群体),则需要将PCA或structure分析的结果作为后续关联分析的协变量,校正它们对关联分析的影响。 运用3:推断进化关系例如下图这篇葡萄群体研究的文章,研究的葡萄品种来源三个地域。绿色的西部葡萄和红色的东部葡萄区分比较明显,而蓝色的中部葡萄夹杂在东、西两个亚群间,和两个亚群有大量重叠。作者从中推断,东、西两个地域的葡萄都有传播到中部地区,并伴随大量杂交,导致中部地区的品种系谱比较混杂,并没有形成自己独立的亚群。 葡萄亚群体的基因混杂现象 PCA分析实操 前期准备 给标记加上IDSNP data通常都是以VCF格式文件呈现,拿到VCF文件的第一件事情就是添加各个SNP位点的ID。 先看一下最开始生成的VCF文件: 原始VCF文件可以看到,ID列都是".",需要我们自己加上去。我用的是某不知名大神写好的perl脚本,可以去我的github上下载,用法: perl path2file/VCF_add_id.pl YourDataName.vcf YourDataName-id.vcf`当然也可以用excel手工添加。添加后的文件如下图所示(格式:CHROMID__POS): 添加ID后VCF文件 SNP位点过滤(Missing rate and maf filtering)SNP位点过滤前需要问自己一个问题,我的数据需要过滤吗? 一般要看后期是否做关联分析(GWAS);如果只是单纯研究群体结构建议不过滤,因为过滤掉低频位点可能会改变某些样本之间的关系;如果需要和表型联系其来做关联分析,那么建议过滤,因为在后期分析中低频位点是不在考虑范围内的,需要保持前后一致。 如果过滤,此处用到强大的plink软件,用法: plink --vcf YourDataName-id.vcf --maf 0.05 --geno 0.2 --recode vcf-iid -out YourDataName-id-maf0.05 --allow-extra-chr参数解释:--maf 0.05:过滤掉次等位基因频率低于0.05的位点;--geno 0.2:过滤掉有20%的样品缺失的SNP位点;--allow-extra-chr:我的参考数据是Contig级别的,个数比常见分析所用的染色体多太多,所以需要加上此参数。 格式转换将vcf文件转换为bed格式文件。 这里注意一点!!!!:应该是软件的问题,需要把染色体/contig名称变成连续的数字(1 to n),不然会报错无法算出结果!(坑) plink --vcf YourDataName-id-maf0.05.vcf --make-bed --out snp --chr-set 29 no-xy参数解释:--chr-set 给出染色体/contig的数目;no-xy 没有xy染色体。 用gcta做PCA分析 gcta输出grm阵列(genetic relationship matrix) gcta64 --make-grm --out snp.gcta --bfile snp --autosome-num 29参数解释:--autosome-num常染色体数目。 gcta计算PCA gcta64 --grm snp.gcta --pca 20 --out snp.gcta参数解读:--pca 20 保留前20个PCA。 特征值结果储存在snp.gcta.eigenval中,特征向量储存在snp.gcta.eigenvec中。 结果处理将特征值结果和特征向量结果用R处理为可读性结果。写好的R包我放在了Github中:PCA2normal_format.R,大家自行下载使用。 如果不想下载,直接复制如下代码: eigvec |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |