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 輸出結(jié)果保留頭文件
#--out pca 輸出結(jié)果文件的前綴
#--allow-extra-chr 染色體設(shè)置
#結(jié)果文件為pca.eigenval,pca.eigenvec

#R畫圖
library(ggplot2)###加載R包ggplot2
library(readr)
library(dplyr)
pca <- read_delim("E:/SNP數(shù)據(jù)處理/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")###讀取各主成分數(shù)據(jù)
data2<-read.table(file="E:/SNP數(shù)據(jù)處理/pca.eigenval",header=F,sep="\t")###讀取各主成分貢獻度數(shù)據(jù)
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="")（###設(shè)置X軸坐標標簽）
ylab<-paste("PC2(",PC2contri,"%)",sep="")###設(shè)置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