差异基因聚类分析：基因型数据绘制PCA图和聚类分析图

差异基因聚类分析：基因型数据绘制PCA图和聚类分析图基因型数据：图片解释，将每个品种用不同的颜色表示，X坐标是PC1，解释24.9%的变异，Y坐标是PC2，解释10.61%的变异，Z坐标是PC3，解释1.02%的变异。可以看到，三个品种在PCA图里面分的比较开，C品种的有两个A和B的点，应该是异常数据。2-D PCA图：图片解释，将每个品种用不同的颜色表示，同时绘制置信区间圆圈，X坐标是PC1，解释24.9%的变异，Y坐标是PC2，解释10.61%的变异。可以看到，三个品种在PCA图里面分的比较开，C品种的有两个A和B的点，应该是异常数据。3-D PCA图：

用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.ped412 re2.ped

SNP个数为：41013

$ wc -l re2.map41013 re2.map

计算思路：

1，对数据进行清洗，将其转化为0 1 2的形式

2，计算G矩阵

3，计算PCA的特征向量和特征值

4，根据特征值计算解释百分比

5，根据特征向量和品种标签，进行PCA的绘制

绘制代码如下：

首先，使用plink命令，将基因型数据转化为012的raw格式：

plink --file re2 --recodeA

结果生成plink.raw文件。

然后使用R语言，计算PCA，并绘制PCA图。

读取数据m012 = fread("plink.raw")
# 保留FID，IID和基因型数据g012 = m012[ -c(3:6)]dim(g012)fid = g012$FIDiid = g012$IIDlibrary(sommer)
# 整理格式，计算G矩阵setDF(g012)rownames(g012) = g012$IIDg012$IID = g012$FID = 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个，用于聚类

n = sample(1:384 100)Gmat1 = Gmat[n n]library(ggtree)library(amap)clu <- hclusterpar(Gmat1)

默认的作图，太挤。

ggtree(clu) geom_tiplab2(offset=0.3)

画圈图：

ggtree(clu layout = "circular") geom_tiplab2(offset=0.1)

可以看到，A，B，C三个品种，可以分开，分得很开。

如果进行GWAS分析，PCA加进去就很有必要！