plink --vcf test.vcf --pca 47 header tabs var-wts --out pca --allow-extra-chr
#--vcf test.vcf 群体SNP的VCF文件
#--pca 47 进行PCA分析，保留前47个主成分
#header tabs var-wts 输出结果保留头文件
#--out pca 输出结果文件的前缀
#--allow-extra-chr 染色体设置
#结果文件为pca.eigenval,pca.eigenvec
#R画图
library(ggplot2)###加载R包ggplot2
library(readr)
library(dplyr)
pca <- read_delim("E:/SNP数据处理/pca.eigenvec",
"\t", escape_double = FALSE, trim_ws = TRUE)
pca <-dplyr::select(pca,FID,PC1,PC2,PC3,PC4)###提取这几列
PCA_DATA <- left_join(pca,group,by = c('FID'='code'))###group文件是分组文件
data<-read.table(file="PCA_DATA",header=T,sep="\t")###读取各主成分数据
data2<-read.table(file="E:/SNP数据处理/pca.eigenval",header=F,sep="\t")###读取各主成分贡献度数据
data3<-data2[c(1:(nrow(data2)-1)),]
PC1contri<-round(data3[1]*100/sum(data3),digits=2)###PC1的贡献率
PC2contri<-round(data3[2]*100/sum(data3),digits=2)###PC2的贡献率xlab<-paste("PC1(",PC1contri,"%)",sep="")（###设置X轴坐标标签）
ylab<-paste("PC2(",PC2contri,"%)",sep="")###设置Y轴坐标标签
p <- ggplot(PCA_DATA,aes(x=PC1,y=PC2,color=number_of_ear_rows,shape=naked_or_hulled)) + geom_point() + labs(x=xlab,y=ylab)
p
##将naked改成hulless
colnames(PCA_DATA) <- c("IID","PC1","PC2","PC3","PC4","number_of_ear_rows","hulless_or_hulled")#改列名
PCA_DATA$hulless_or_hulled[which(PCA_DATA$hulless_or_hulled =='naked')]<- 'hulless'#把hulless or hulled 那一列的naked改成hulless
p <- ggplot(PCA_DATA,aes(x=PC1,y=PC2,color=number_of_ear_rows,shape=hulless_or_hulled)) + geom_point() + labs(x=xlab,y=ylab)#画图
p
