基因型数据的聚类分析和PCA分析

科技   科技   2024-10-08 20:05   河南  

大家好,我是邓飞。


基因型数据的聚类分析和主成分分析(PCA)是常用的数据降维和探索性数据分析方法。


1,聚类分析

目标:聚类分析旨在将样本分成若干组,使得组内样本的相似度最大,组间样本的相似度最小。聚类分析可以帮助识别群体结构、发现自然分组、识别相似基因型的个体群体等。
方法:
  • K-均值聚类(K-means):根据预设的簇数(K),迭代地分配样本点到最近的簇中心,然后更新簇中心的位置。
  • 层次聚类(Hierarchical Clustering):创建一棵树状的层次结构,样本逐步合并成簇或从簇中拆分。

2,主成分分析(PCA)

目标:PCA用于将高维数据降维到低维空间,同时尽量保留原始数据的方差信息。

方法:
  • 特征值分解或奇异值分解:计算协方差矩阵并进行特征值分解,从中选择最大的特征值对应的特征向量作为主成分。

步骤:
  • 标准化数据(中心化和缩放)。
  • 计算协方差矩阵。
  • 计算协方差矩阵的特征值和特征向量。
  • 根据特征值排序选择前几个特征向量作为主成分。

应用:
PCA在遗传学中常用于可视化高维基因型数据、识别样本间的变异模式、纠正群体分层等。


3,GWAS中的主成分分析(PCA)


用PCA作为GWAS的协变量,相当于将品种结构考虑进去。它类似将不同品种作为协变量,或者将群体结构矩阵Q作为协变量。


很多软件可以分析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 = NULLg012$FID = NULLGmat = 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"))


4,基因型数据的聚类分析


聚类分析思路:


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三个品种,可以分开,分得很开。

5,PCA分析和聚类分析的区别

聚类分析:结果解释集中在样本的分组和聚类特性上,适合用于定性分析样本之间的相似性和差异性。

PCA分析:结果解释集中在主成分的线性组合及其能够解释的总方差上,适合用于理解数据的主要变异和趋势,以及对噪声进行过滤。

从计算方法上来看,两者都是先计算亲缘关系矩阵,然后对亲缘关系G矩阵,进行进一步的分析而来。

推荐阅读:

想要更好的学习和交流,快来加入飞哥的知识星球,这是一个生物统计+数量遗传学+GWAS+GS的社区,在这里你可以向飞哥提问、帮你制定学习计划、跟着飞哥一起做实战项目,冲冲冲。点击这里加入吧:飞哥的学习圈子


1,快来领取 | 飞哥的GWAS分析教程


2,飞哥汇总 | 入门数据分析资源推荐


3,数量遗传学,分享几本书的电子版


4,R语言学习看最新版的电子书不香嘛?


5,书籍及配套代码领取--统计遗传分析导论


6,十一在家把GWAS分析学会吧!


育种数据分析之放飞自我
本公众号主要介绍动植物育种数据分析中的相关问题, 算法及程序代码.
 最新文章