群体遗传专题-PCA分析
今周开始给大家逐步介绍一些群体遗传分析会用到的生物信息学分析方法。所介绍的大部分分析都是基于你已经知道怎样call snp，然后基于snp的结果，进行一系列群体遗传相关的分析。
对于那些对SNP变异还不熟悉的朋友，在阅读本文之前，我给你推荐一个由咱们生信技能树资深大神写的博文：我见过整理的最好最易懂的最完整的SNP calling的博文。
PCA概念介绍
PCA 全称是Principal Component Analysis，又叫做主成分析。是一种分析、简化数据集的技术。主成分分析经常用于减少数据集的维数，同时保持数据集中的对方差贡献最大的特征。这是通过保留低阶主成分，忽略高阶主成分做到的。这样低阶成分往往能够保留住数据的最重要方面。
主成分分析由卡尔·皮尔逊于1901年发明，用于分析数据及建立数理模型。其方法主要是通过对协方差矩阵进行特征分解，以得出数据的主成分（即特征向量）与它们的权值。PCA是最简单的以特征量分析多元统计分布的方法。其结果可以理解为对原数据中的方差做出解释：哪一个方向上的数据值对方差的影响最大？换而言之，PCA提供了一种降低数据维度的有效办法。
PCA基本原则在最大程度反映原变量所代表的信息，同时保证新变量之间信息不重复。在生物学上，常常用于将SNP信息浓缩为几个新的变量。
PCA应用
群体分层分析和推断进化关系
可以准确的展示群体分层的类型，一般可以与phylogentic tree，structure，admixture的结果互相验证。
如图左边绿色是野生的大麦，右边橙色是驯化的大麦，中间蓝色是admixture混合的大麦。不同品系的大麦可以看到明显的分层，并且具有野生和驯化大麦特征的混合小麦位于野生还有驯化小麦之间。很好的反映了，混合小麦是其两者的一种中间过度形态，具有两者的特征。
2.检查outliner
在一个群体中，如果其中有一些样本取样错误，或者测序时候有严重的污染，这样往往会导致outliner的产生，如下图，右上方的那个样本就是明显的outlier，有着与其他样本不同，偏高的PCA比值。当确认outliner后，可以准确校正群体的样本，减少对后续关联分析的影响。
实战
这里会用到一组GBS的数据，由于本节内容主要讲述PCA，所以就不详细讲解具体数据结构还有内容（下章再回来讲）。
Load 对应的R包
library(vcfR)
library(poppr)
library(ape)
library(RColorBrewer)
读取vcf文件，还有数据信息文件，将vcf转化成genlight格式
rubi.VCF <- read.vcfR("prubi_gbs.vcf")
rubi.VCF
pop.data <- read.table("population_data.gbs.txt", sep = "\t", header = TRUE)
pop(gl.rubi) <- pop.data$State
gl.rubi <- vcfR2genlight(rubi.VCF)
使用glPCA function进行PCA分析
rubi.pca <- glPca(gl.rubi, nf = 3)
分析完后，对特征值用图形进行大概了解。
variance <- as.matrix(rubi.pca$eig)
par(mar = c(5.1,4.5,4.1,4.5))
plot(variance[1:10], type = "b", col = "blue", pch = 20, xaxt="n",yaxt="n", xlab="PCs",ylab="Variance")
axis(side = 1, at = 1:10)
axis(side = 2, at = round(seq(0, max(variance), length = 12)))
ec <- vector(mode="numeric", length=length(variance))
for (i in 1:length(variance)){
ec[i] <- sum(variance[1:i])/sum(variance)
}
par(new=TRUE)
plot(ec[1:10], type="b", col="red",xaxt="n",yaxt="n",xlab="",ylab="", ylim = c(0,1), pch = 20)
axis(side = 4, at = seq(0, 1, by = 0.1))
mtext("Cumulative variance explained", side=4, line=3)
legend(x = 5, y = 1.0, bty = "n",col=c("blue", "red"),lty=1,legend=c("Variance","Cumulative variance explained"))
结果如下图：前3个PCA累计解释63.183％的数据差异。
要查看PCA的结果，我们可以使用R包ggplot2。我们需要将包含主要组件（rubi.pca $scores）的数据框转换为新的对象rubi.pca.scores。另外，我们将在rubi.pca.scores对象中添加人口值作为新列，以便能够按种群对样本着色。ggplot2将绘制PCA图，根据总体对样本着色，并创建包含每个群体95％数据的椭圆圈。
rubi.pca.scores <- as.data.frame(rubi.pca$scores)
pop(gl.rubi) <- pop.data$State
rubi.pca.scores$pop <- pop(gl.rubi)
cols <- brewer.pal(n = nPop(gl.rubi), name = "Dark2")
set.seed(9)
p <- ggplot(rubi.pca.scores, aes(x=PC1, y=PC2, colour=pop))
p <- p + geom_point(size=2)
p <- p + stat_ellipse(level = 0.95, size = 1)
p <- p + scale_color_manual(values = cols)
p <- p + geom_hline(yintercept = 0)
p <- p + geom_vline(xintercept = 0)
p <- p + theme_bw()
p
结果如下：PCA分析将WA和CA样本（左边）还有样本CA和WA（右边）分离开来。然后样本CA比较显著的聚集在一起。
本文的教程还有数据都是基于这个网站上的内容：
http://grunwaldlab.github.io/Population_Genetics_in_R/index.html
有条件有兴趣的朋友也可以自行进行更多的深入阅读。
猜你喜欢
生信菜鸟团-专题学习目录（6）
还有更多文章，请移步公众号阅读
▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。
▼ 如果你是初学者，请关注下面的生信菜鸟团，了解生信基础名词，概念，扎实的打好基础，争取早日入门。