library(SNPRelate) library(ggplot2) vcf.fn <- "spipsforPCA.recode.vcf" snpgdsVCF2GDS(vcf.fn, "all.recode.SPIPS.gds", method="biallelic.only") snpgdsSummary("all.recode.SPIPS.gds") genofilefilter_allrecode <- snpgdsOpen("all.recode.SPIPS.gds") pcafilter_allrecode <- snpgdsPCA(genofilefilter_allrecode, autosome.only=FALSE, num.thread=2) write.csv(pcafilter_allrecode$eigenval,"eigen.csv") tabfilter_allrecode <- data.frame(sample.id = pcafilter_allrecode$sample.id, EV1 = pcafilter_allrecode$eigenvect[,1], # the first eigenvector EV2 = pcafilter_allrecode$eigenvect[,2], # the second eigenvector stringsAsFactors = FALSE) plot(tabfilter_allrecode$EV1, tabfilter_allrecode$EV2, xlab="eigenvector 1", ylab="eigenvector 2") write.csv(tabfilter2_allrecode$sample.id,file = 'samplenames_SPIPS.csv') write.csv(tabfilter_allrecode$EV1,file="PC1.csv") sample.id <-pcafilter_allrecode$sample.id sample.id # write.csv(sample.id, file = 'sampleid.csv') pop_code <- scan("popSPIPS.group", what=character()) write.csv(cbind(sample.id, pop_code),"~/Desktop/names.csv") tabfilter2_allrecode <- data.frame(sample.id = pcafilter_allrecode$sample.id, pop = factor(pop_code)[match(pcafilter_allrecode$sample.id, sample.id)], EV1 = pcafilter_allrecode$eigenvect[,1], # the first eigenvector EV2 = pcafilter_allrecode$eigenvect[,2], # the second eigenvector stringsAsFactors = FALSE) tabfilter2_allrecode palette(c("Red","Green","Blue")) ggplot(data = tabfilter2_allrecode,aes(x=EV1,y=EV2,colour=pop))+ geom_point(shape=16) + theme_bw() + ggtitle("Single Spaced Plants PCA") + theme(plot.title = element_text(hjust = 0.5))+ labs(x="PC1 (1.22%)",y="PC2 (5.23%)")+ guides(colour=guide_legend(title="Ecotype")) + theme(axis.text=element_text(size=18), axis.title=element_text(size=18,face="bold"), legend.text = element_text(size=18), legend.title=element_text(size=18,face="bold"), plot.title=element_text(size=24,face='bold')) p<-ggplot(tabfilter2_allrecode,aes(x = EV1, y = EV2, colour = pop)) + geom_point() + scale_color_manual(name="Ecotypes", labels=c("Dry Ecotype","Mesic Ecotype","Wet Ecotype"), values=c("#F8766D","#00BA38","#619CFF")) + theme_bw() + labs(x="PC1 (5.23%)", y="PC2 (1.22%)") + theme(legend.position='bottom') + ylim(-0.2,0.2) + xlim(-0.2,0.2) + coord_fixed() p ggsave("PCAGenetic.jpeg", dpi=1200)