根据vcf文件做主成分分析、构建进化树小例子

参考文章

https://zhuanlan.zhihu.com/p/45950835

使用的数据

https://grunwaldlab.github.io/Population_Genetics_in_R/gbs_analysis.html
下载这篇文章中用到的数据

主成分分析代码

library(gdsfmt)
library(SNPRelate)
setwd("../../Population_Genetic_tutorial/Practice_data/")
vcffile<-"prubi_gbs.vcf"
snpgdsVCF2GDS(vcffile,"new_geno.gds",method="biallelic.only")
genofile<-snpgdsOpen("new_geno.gds")
pca<-snpgdsPCA(genofile,autosome.only = F)
pc.percent<-pca$varprop*100
pcadf<-data.frame(sampleid=pca$sample.id,
                  PC1=pca$eigenvect[,1],
                  PC2=pca$eigenvect[,2],
                  stringsAsFactors = F)
plot(pcadf$PC1,pcadf$PC2)
head(pcadf)
popdf<-read.table("population_data.gbs.txt",
                  sep="\t",header=T)
head(popdf)
pcadata<-cbind(pcadf,popdf)
head(pcadata)
which(pcadata$sampleid != pcadata$AccessID)
library(ggplot2)
ggplot(pcadata,aes(x=PC1,y=PC2,color=State))+
  geom_point(aes(shape=State))+
  theme_bw()

image.png

聚类分析代码

ibs.hc<-snpgdsHCluster(snpgdsIBS(genofile,num.thread=2,
                                 autosome.only = F))
rv<-snpgdsCutTree(ibs.hc)
pdf(file="tree.pdf",width=15)
plot(rv$dendrogram,leaflab="perpendicular")
dev.off()

tree.png

如何美化这个结果还得好好学习！
欢迎大家关注我的公众号
小明的数据分析笔记本

公众号二维码.jpg