用PCA做为GWAS的协变量，相当于将品种结构考虑进去。它类似将不同品种作为协变量，或者将群体结构矩阵Q作为协变量。
下面看一下利用基因型SNP数据进行PCA计算，以及可视化的分析。
很多软件可以分析PCA，这里介绍一下使用plink软件和R语言，进行PCA分析，并且使用ggplot2绘制2D和3D的PCA图。
绘制后的图如下：
2-D PCA图：
图片解释，将每个品种用不同的颜色表示，同时绘制置信区间圆圈，X坐标是PC1，解释24.9%的变异，Y坐标是PC2，解释10.61%的变异。可以看到，三个品种在PCA图里面分的比较开，C品种的有两个A和B的点，应该是异常数据。
3-D PCA图
图片解释，将每个品种用不同的颜色表示，X坐标是PC1，解释24.9%的变异，Y坐标是PC2，解释10.61%的变异，Z坐标是PC3，解释1.02%的变异。可以看到，三个品种在PCA图里面分的比较开，C品种的有两个A和B的点，应该是异常数据。
基因型数据：
共有3个品种A，B，C，共有412个个体。其中：
- A品种有200个体
- B品种有100个体
- C品种有112个体
$ wc -l re2.ped
412 re2.pedSNP个数为：41013
代码语言：javascript复制$ wc -l re2.map
41013 re2.map计算思路：
1，对数据进行清洗，将其转化为0,1,2的形式
2，计算G矩阵
3，计算PCA的特征向量和特征值
4，根据特征值计算解释百分比
5，根据特征向量和品种标签，进行PCA的绘制
绘制代码如下：
首先，使用plink命令，将基因型数据转化为012的raw格式：
代码语言：javascript复制plink --file re2 --recodeA结果生成plink.raw文件。
然后使用R语言，计算PCA，并绘制PCA图。
代码语言：javascript复制读取数据
m012 = fread("plink.raw")
# 保留FID，IID和基因型数据
g012 = m012[,-c(3:6)]
dim(g012)
fid = g012$FID
iid = g012$IID
library(sommer)
# 整理格式，计算G矩阵
setDF(g012)
rownames(g012) = g012$IID
g012$IID = NULL
g012$FID = NULL
Gmat = A.mat(g012-1)
# 计算特征值和特征向量
re = eigen(Gmat)
# 计算解释百分比
por = re$values/sum(re$values)
# 整理格式
pca_re1 = re$vectors[,1:3]
pca_re2 = data.frame(pca_re1,Ind = iid)
pca_re2$Gen = fid
# 把PCA1，PC2，PC2的j解释百分比，命名为相应的轴
xlab = paste0("PC1(",round(por[1]*100,2),"%)")
ylab = paste0("PC2(",round(por[2]*100,2),"%)")
zlab = paste0("PC3(",round(por[3]*100,2),"%)")
# 绘制2-D PCA图
ggplot(pca_re2, aes(x=X1, y=X2,color=Gen)) geom_point(size=2)
# stat_ellipse(level = 0.95, size = 1)
stat_ellipse(aes(fill=Gen),
type ="norm", geom ="polygon",alpha=0.2,color=NA)
geom_hline(yintercept = 0) # 添加x坐标
geom_vline(xintercept = 0) # 添加y坐标
labs(x = xlab,y = ylab,color="")
guides(fill=F)
theme_bw() # 主题
# 绘制3-D PCA图
library(scatterplot3d)
pca_re2 = pca_re2 %>% mutate(colour = case_when(
Gen == "A" ~ "red",
Gen == "B" ~ "green",
Gen == "C" ~ "blue",
))
scatterplot3d(pca_re2[,1:3],color=pca_re2$colour,
pch = 16,angle=30,
box=T,type="p",
xlab = xlab, ylab = ylab, zlab = zlab,main = "3D PCA Plot",
lty.hide=2,lty.grid = 2)
legend("topright",c("A","B","C"),fill=c('red','green',"blue"))
聚类分析思路：
1，计算个体的亲缘关系矩阵G矩阵
2，因为数据太多，随机选择100个，用于聚类
代码语言：javascript复制n = sample(1:384,100)
Gmat1 = Gmat[n,n]
library(ggtree)
library(amap)
clu <- hclusterpar(Gmat1)默认的作图：
代码语言：javascript复制ggtree(clu) geom_tiplab2(offset=0.3)
画圈图：
代码语言：javascript复制ggtree(clu,layout = "circular") geom_tiplab2(offset=0.1)
可以看到，A，B，C三个品种，可以分开，分得很开。
